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Abstract 

This  thesis  proposes  a  methodology  for  optimally  sampling  a  chemical  hazard 
area  subsequent  to  a  chemical  weapons  attack.  The  objective  is  to  identify  the  max¬ 
imum  number  of  areas  that  no  longer  require  protective  gear  for  safe  operations.  We 
model  the  area  as  an  undirected  graph  and  employ  network  analysis  techniques  to 
provide  a  methodological  framework  for  identifying  an  optimal  sampling  sequence 
within  a  fixed  time  limit.  We  propose  four  models  that  characterize  the  secondary 
vapor  concentrations:  i)  static  and  deterministic,  ii)  static  and  stochastic,  iii)  dy¬ 
namic  and  deterministic,  and  iv)  dynamic  and  stochastic.  Comparisons  of  the  static 
cases  and  their  dynamic  counterparts  demonstrate  the  impact  of  temporal  evolution 
of  vapor  concentrations  on  the  optimal  sampling  path.  We  conclude  that  the  num¬ 
ber  of  safe  areas  may  be  either  under-  or  over-estimated  depending  on  the  assumed 
nature  of  the  secondary  vapors. 
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OPTIMAL  SAMPLING  OF  A  CHEMICAL 
HAZARD  AREA 


1.  Introduction 

1.1  Background 

The  willingness  of  contemporary  terrorist  organizations  to  use  asymmetric  mea¬ 
sures  against  their  adversaries  suggests  the  need  for  defending  U.S.  military  installa¬ 
tions  and  interests  against  chemical  or  biological  attacks.  According  to  [6] ,  chemical 
and  biological  weapons  (CBWs)  are  considered  to  be  weapons  of  mass  destruction 
(WMD)  because  they  are  capable  of  “high  order  destruction  and/or  of  being  used  in 
such  manner  as  to  destroy  large  numbers  of  people.”  The  devastation  CBWs  may 
cause,  and  the  ease  with  which  they  may  now  be  developed,  stored,  and  transported 
has  made  protecting  against  them  a  primary  concern  for  military  installations.  In  the 
event  of  an  attack  on  any  United  States  military  installation  or  interest,  a  method 
for  identifying  the  most  hazardous  areas  is  required  to  protect  human  life  and  to 
restore  operations  to  full  capability  as  soon  as  possible. 

The  terrorist  attacks  on  the  United  States  on  September  11,  2001  have  dra¬ 
matically  altered  the  views  of  the  American  public  and  militaries.  As  noted  in  [6], 
the  threats  to  U.S.  targets,  military  combat  personnel,  and  American  military  in¬ 
stallations  in  foreign  countries  are  very  credible.  However,  steps  have  been  taken 
to  reduce  the  threat  of  chemical  or  biological  attacks  on  groups  or  nations.  The 
Chemical  Weapons  Convention  (CWC)  in  1997,  and  Biological  Weapons  Convention 
(BWC)  were  “adopted  to  stifle  proliferation  of  chemical  and  biological  weapons.” 
However,  there  are  shortcomings  with  the  treaties.  They  do  not  focus  on  the  small 
quantities  that  a  terrorist  group  may  employ,  but  on  the  large  developmental  pro- 
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grams  of  entire  nations.  Also,  they  do  not  regulate  “dual  use”  items  or  chemicals 
such  as  Toxic  Industrial  Chemicals  (TICs)  that  can  double  as  weapons.  These  in¬ 
clude  chemicals  such  as  methyl  isocynate  which  was  used  in  Bhopal,  India  (Dec  1984) 
that  are  not  included  in  the  schedules  of  chemicals  listed  in  the  CWC  [6].  Due  to 
these  shortcomings,  emergency  response  and  assessment  teams  must  be  prepared  to 
act  following  an  attack. 

The  Department  of  Defense  Chemical,  Biological,  Radiological,  and  Nuclear 
Defense  Program  (CBRNDP)  Annual  Report  to  Congress  dated  May  2004  [5]  con¬ 
tains  significant  literature  on  the  continuing  need  to  defend  against  weapons  of  mass 
destruction  (WMD)  today.  Among  several  other  topics,  the  report  includes  infor¬ 
mation  from  an  unclassified  CIA  report  to  Congress  concerning  the  acquisition  of 
technology  that  relates  to  WMD  and  advanced  conventional  munitions  which  states 
current  threats  to  the  United  States  and  its  assets.  It  also  focuses  on  contamination 
avoidance  as  well  as  improving  the  individual  protective  equipment  for  the  personnel 
who  may  be  exposed  to  contamination. 

According  to  the  CBRNDP  report,  the  CIA  has  provided  evidence  that  several 
countries  have  an  active  interest  in  WMD.  These  countries  include  Iran,  Iraq,  North 
Korea,  Libya,  Syria,  Sudan,  India,  and  Pakistan  and  some  of  the  key  supplier  coun¬ 
tries  are  reported  to  be  Russia,  China,  and  North  Korea.  Disturbingly,  many  of  the 
33  “designated  terrorist  organizations”  and  other  non-state  actors  across  the  world 
have  shown  interest  in  CBRN  [5].  With  all  of  these  countries  and  organizations  ex¬ 
pressing  interest  in  WMD,  it  is  prudent  for  the  LhS.  to  be  concerned  with  defending 
against  these  weapons  and  to  develop  strategic  policies  that  can  be  implemented  in 
the  event  of  such  an  attack  to  protect  human  life  and  to  restore  hazardous  areas  to 
full  operational  capabilities  as  quickly  as  possible. 

Subsequent  to  a  successful  chemical  attack  on  any  fixed  operational  site,  miliary 
or  otherwise,  vapor  contamination  can  hinder  operations  for  possibly  a  significant 
amount  of  time.  The  contamination  results  from  persistent  chemical  agents  that 
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deposit  on  various  surfaces  and  then  evaporate  into  the  air,  creating  a  hazardous 
vapor  cloud.  Atmospheric  conditions  may  then  cause  this  cloud  to  travel  to  parts  of 
the  site  that  may  not  have  been  directly  hit  and  cause  contamination,  threatening  the 
health  of  personnel  in  the  area  as  well  as  their  productivity.  Decision  makers  must 
determine  the  level  of  mission-oriented  protective  posture  (MOPP)  for  personnel  in 
various  work  areas.  The  protective  gear  donned  by  personnel  in  (potentially)  infected 
areas  can  be  cumbersome  and  cause  heat  stress;  therefore,  commanders  must  take 
into  account  the  specific  threat,  environmental  conditions  that  affect  the  evaporation 
rates,  work  rates,  the  level  of  task  difficulty,  and  mission  requirements  when  deciding 
the  appropriate  MOPP  level  at  any  point  in  time  following  a  chemical  attack. 

It  is  imperative  to  identify  the  areas  of  the  site  which  can  safely  return  to  full 
operational  capabilities  as  soon  as  possible.  Knowing  which  areas  of  the  site  crews 
should  sample  first  in  order  to  identify  where  gear  is  (or  is  not)  required  is  essential 
to  restoring  the  site  to  full  capabilities.  At  each  sampled  area,  instrument  readings 
must  be  obtained  so  that  secondary  vapor  concentrations  may  be  assessed. 

Developing  a  strategy  for  prioritizing  the  areas  to  be  searched  that  identifies 
the  maximum  number  of  areas  in  their  optimal  level  of  safety  is  the  focus  of  this 
thesis.  That  is,  by  carefully  choosing  which  areas  to  search  to  obtain  actual  vapor 
readings,  a  maximum  number  of  areas  can  be  reached  within  an  allotted  time  while 
only  searching  those  that  are  “critical.”  Critical  areas  are  identified  as  those  where 
the  vapor  concentration  level  has  decreased  below  a  fixed  threshold  that  dictates  the 
need  for  protective  gear. 

The  basic  problem  can  be  viewed  as  a  set  covering  problem.  There  exists  a  set 
of  areas  that  need  to  be  reached  from  some  starting  location  while  minimizing  (or 
maximizing)  some  quantity.  In  this  case,  we  seek  to  cover  the  maximum  number  of 
areas  where  protective  gear  is  not  required.  This  problem  increases  in  complexity 
due  to  the  addition  of  constraints,  specifically  a  time  constraint,  forcing  a  feasible 
solution  to  the  real  world  problem  to  satisfy  these  time  restrictions.  For  example, 
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the  solution  must  be  a  path  through  the  site  that  does  not  allow  subtours  (disjoint 
paths)  where  the  crew  cannot  move  from  one  path  to  another. 

The  approach  used  in  this  thesis  has  not  been  previously  applied  to  a  problem 
considering  the  identification  of  secondary  vapor  contamination  to  an  operational 
site  subsequent  to  a  chemical  weapon  attack.  We  model  this  problem  as  a  network 
and  employ  operations  research  techniques  to  solve  it.  This  approach  also  addresses 
the  spatial  and  temporal  dynamics  of  the  real  problem  by  incorporating  them  into 
the  model.  Though  extensive  research  has  been  ongoing  for  defending  against  a 
chemical  weapons  attack,  we  focus  our  attention  on  obtaining  a  policy  for  the  full 
restoration  of  operational  capabilities  as  soon  as  possible  should  a  chemical  attack 
occur. 

1.2  Problem  Definition  and  Methodology 

This  thesis  is  primarily  concerned  with  the  development  of  a  methodology  for 
incorporating  available  vapor  level  data  to  “intelligently”  assess  the  actual  threat  to 
personnel  in  areas  throughout  a  fixed  military  installation.  Thus,  the  problem  of 
optimally  sampling  a  chemical  hazard  area  is  our  focus. 

The  dynamics  of  the  problem  may  be  described  as  follows.  Following  a  chemical 
attack,  detectors  such  as  M-8  or  M-9  paper  are  available  for  reading.  They  contribute 
a  “yes”  or  “no”  response  as  to  whether  or  not  a  liquid  drop  of  chemical  has  touched 
the  paper.  As  the  liquid  begins  to  evaporate  and  enter  the  atmosphere,  the  wind 
carries  it  downstream  and  begins  to  disperse  the  hazardous  vapors  to  other  areas. 
Additionally,  as  time  progresses,  the  concentration  of  the  vapors  begins  to  diminish 
and  some  areas  may  become  safer.  It  is  important  to  know  when  the  areas  begin 
to  return  to  normal  so  that  personnel  may  return  to  normal  operations.  Thus, 
methods  for  searching  the  site,  with  the  objective  of  identifying  the  areas  with  a 
vapor  concentration  below  a  specific  threshold,  are  considered. 
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This  problem  displays  characteristics  of  several  well-known  problems.  Part  of 
the  objective  is  the  same  as  for  the  Travelling  Salesperson  Problem  (TSP),  where  we 
seek  to  minimize  the  time  to  reach  all  cities  on  a  tonr.  ffowever,  for  this  problem, 
the  “salesperson”  has  a  fixed  amount  of  time  to  reach  all  of  his  or  her  destinations. 
Thus,  it  is  possible  that  not  all  “cities”  will  be  reached.  Also  we  seek  to  maximize 
the  reward  for  visiting  “cities”  and  the  “salesperson”  does  not  have  to  return  to  the 
starting  location. 

Another  well-known  network  problem  similar  to  the  one  addressed  here  is  the 
shortest  path  problem  (SPP).  The  entity  in  the  shortest  path  problem  must  find  the 
path  of  minimum  length  (or  time)  from  a  starting  location  (source)  to  a  specified 
destination  (sink).  Our  problem  strays  from  the  SPP  because  we  do  not  have  a 
distinct  sink  node  and  we  impose  an  upper  bound  on  the  length  (time)  of  any  path 
in  our  problem.  It  also  varies  from  the  SPP  because  we  seek  to  sample  as  many  areas 
as  possible  rather  than  simply  finding  a  path  of  shortest  length.  There  is  a  reward 
associated  with  reaching  a  node  in  our  problem  which  is  not  necessarily  included  in 
the  standard  SPP. 

This  problem  can  also  be  viewed  as  a  knapsack  problem  where  a  fixed  amount 
of  resource  is  available  for  use  in  maximizing  the  objective.  In  our  problem,  the 
resource  is  time,  and  we  seek  to  maximize  the  reward.  Contrary  to  the  knapsack 
problem,  the  order  in  which  the  items  are  placed  in  the  “knapsack”  is  important. 

Since  we  seek  to  reach  all  the  nodes,  if  possible,  this  problem  shares  a  similarity 
with  the  minimum  spanning  tree  problem.  However,  in  the  minimum  spanning  tree, 
each  node  must  be  connected  but  they  may  not  form  a  connected  path  through 
the  network.  We  seek  a  path  that  can  be  followed  from  beginning  to  end  without 
reaching  a  dead  end  or  becoming  stuck  in  a  cycle. 

It  is  clear  the  problem  of  this  thesis  shares  characteristics  with  several  well- 
known  optimization  problems.  However,  new  methodologies  must  be  developed  to 
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solve  this  particular  problem.  The  proposed  methodology  incorporates  aspects  of 
each  of  those  problems,  but  will  be  distinct  from  them  in  its  formulation. 

Our  main  contribution  to  the  operations  research  literature  is  an  optimal  sam¬ 
pling  procedure  for  a  single  crew  to  assess  the  current  threat  level  at  a  fixed  opera¬ 
tional  site  in  order  to  reduce  MOPP  levels  as  quickly  as  possible.  We  consider  four 
distinct  cases  of  this  problem.  First,  we  consider  the  case  where  vapor  concentra¬ 
tion  levels  are  static  and  deterministic.  Next,  we  assume  the  vapor  concentrations 
evolve  temporally,  but  remain  deterministic.  Third,  we  consider  stationary  proba¬ 
bility  distributions  for  the  vapor  concentrations  and  finally,  we  consider  the  case  in 
which  vapor  concentrations  evolve  according  to  a  time-variant  probability  distribu¬ 
tion  in  each  area.  The  methodologies  as  stated  are  progressively  more  complex.  The 
ultimate  solution  is  a  time-adaptive  policy  directing  the  search  crew  through  the 
site  to  maximize  the  number  of  areas  identified  to  have  vapor  concentrations  below 
the  threshold.  Thus,  the  maximum  number  of  areas  have  fully  restored  operational 
capabilities. 

1 . 3  Thesis  Outline 

Chapter  2  examines  current  literature  in  the  areas  of  the  ongoing  threat  of 
biological  and  chemical  weapons  to  U.S.  installations,  domestic  and  foreign,  as  well 
as  U.S.  interests,  and  steps  taken  to  defend  against  such  an  attack.  Methodologies 
for  formulating  and  solving  problems,  such  as  the  one  in  this  thesis,  are  also  ex¬ 
amined  and  discussed.  In  chapter  3,  we  formally  define  the  problem  and  provide  a 
mathematical  programming  formulation  for  each  of  the  cases  mentioned  above.  The 
first  section  describes  the  network  model  definitions  and  the  second  section  states 
the  model  assumptions  and  defines  each  of  the  four  cases.  Chapter  4  presents  nu¬ 
merical  illustrations  for  each  of  the  four  cases  and  compares  the  results  of  the  static 
cases  with  their  dynamic  counterparts.  Finally,  Chapter  5  provides  some  concluding 
remarks,  recommendations,  and  directions  for  extensions  and  future  research. 
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2.  Relevant  Literature 


This  chapter  reviews  relevant  contributions  to  the  area  of  emergency  response 
subsequent  to  a  biological  or  chemical  attack,  and  a  review  of  literature  on  network 
problems.  Specifically,  literature  regarding  stochastic,  time-varying  networks,  time- 
adaptive  route  choice  problems,  and  shortest  path  problems  are  reviewed. 

2.1  Some  Historical  Perspectives 

The  use  of  biological  weapons  is  not  a  modern  idea.  It  has  existed  since  600 
B.C.  when  the  Athenians  were  known  to  contaminate  rivers  with  skunk  cabbage  to 
incapacitate  their  enemies  with  a  violent  sickness  [2],  During  the  1300s,  Tartarians 
utilized  corpses  infected  with  the  plague  by  catapulting  them  over  the  city  walls 
of  Kaffa,  possibly  beginning  the  Black  Death  [19]  and  in  America,  blankets  con¬ 
taminated  with  smallpox  were  distributed  to  some  Native  American  tribes  in  the 
1700s  [19].  Most  recently,  the  2001  introduction  of  Anthrax  into  the  United  States 
Postal  Service  reminded  the  American  government  and  domestic  population  just  how 
serious  biological  weapons  are  [20]  and  that  they  are  a  contemporary  threat. 

The  first  known  manufactured  chemical  weapons  date  back  to  1823  when  mus¬ 
tard  gas  (levinstein  mustard)  was  first  synthesized.  The  first  use  of  a  chemical 
weapon  dates  back  to  1914  when  mustard  gas  was  used  during  World  War  I  by  Ger¬ 
man  forces  [10].  More  recently,  the  1995  attack  on  a  Tokyo  subway  by  a  Japanese 
cult  employed  sarin  nerve  gas  killing  twelve  people  and  injuring  thousands  [3].  More 
notably,  in  March  1998,  the  Iraqi  government  used  a  “cocktail”  of  at  least  four  chemi¬ 
cal  weapons  (mustard  gas,  sarin,  tabun,  and  VX)  on  Halabja,  killing  more  than  5,000 
individuals  [7]. 

It  is  evident  that,  throughout  history,  both  chemical  and  biological  weapons 
have  wreaked  havoc  on  nations  and  communities.  Currently,  concerns  about  chem¬ 
ical  weapons  development  and  use  are  elevating  for  governments  worldwide  clue  to 
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scientific  and  technological  advancements  and  the  ease  with  which  they  can  be  ob¬ 
tained  by  adversaries  and  terrorists  [6].  Today,  they  can  be  manufactured  in  the 
same  factories  that  industrial  or  agricultural  chemicals  are  made.  Contemporary 
chemical  weapons  are  easy  to  store  as  a  liquid,  which  can  be  dispersed  in  vapor  form 
(type  of  gas),  or  as  an  aerosol,  (suspension  of  tiny  liquid  or  solid  particles  in  a  gas) 
[15].  Chemical  and  biological  weapons  can  be  created  using  minimal  technology; 
thus  making  them  available  to  any  terrorist  group,  country,  or  state  who  desires 
them.  Advanced  nations  can  produce  more  elaborate  weapons  with  the  use  of  chem¬ 
ical  engineering,  pharmaceutical,  or  biotechnology  industries.  More  importantly,  the 
globalization  of  these  industries,  scientific  engineering,  and  technical  personnel  ex¬ 
changes,  have  made  sharing  information  much  easier  leading  to  even  easier  access  to 
chemical  and  biological  technology  [6]. 

Contemporary  terrorist  groups  are  willing  to  use  asymmetric  measures  to 
accomplish  their  goals,  and  therefore,  threaten  the  use  of  chemical  or  biological 
weapons.  Due  to  the  absence  of  dominating  global  powers  other  than  the  United 
States,  the  threat  is  very  real  and  alarming.  Thus,  chemical  and  biological  defense 
has  received  much  attention,  especially  the  development  of  chemical  sensors  and 
instruments  for  detecting  and  quantifying  chemical  vapors.  Other  areas  of  focus 
include  the  development  of  emergency  response  strategies  for  affected  regions  and 
individuals,  and  restoration  of  a  site  subsequent  to  a  chemical  or  biological  attack. 

2.2  Current  and  Developing  Research  Areas 

The  Department  of  Defense  Chemical,  Biological,  Radiological,  and  Nuclear 
Defense  Program  (CBRNDP)  Annual  Report  to  Congress  from  May  2004  [5]  stresses 
contamination  avoidance  including  reconnaissance,  detection,  and  identification.  If 
a  chemical  agent  release  is  unavoidable,  these  three  areas  are  key  to  ensuring  forces 
in  the  affected  area  assume  the  optimal  level  of  mission  oriented  protective  posture 
(MOPP)  to  sustain  operations  and  restore  military  installations  to  full  operational 
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capabilities  as  soon  as  possible.  Currently,  advanced  technologies  are  being  developed 
for  chemical  and  biological  standoff  detection,  early  warning  detection,  miniaturiza¬ 
tion,  and  interconnectivity.  There  are  other  areas  also  receiving  attention  such  as 
enhancements  in  detection  sensitivity,  interference  rejection,  logistics  supportability, 
and  affordability  [5]. 

The  CBRNDP  report  [5]  states  that  contamination  avoidance  seeks  “to  provide 
a  real-time  capability  to  detect,  identify,  characterize,  quantify,  locate,  and  warn 
against  all  known  or  validated  CBRN  warfare  agent  hazards.”  The  optimization 
of  sensor  technologies  is  being  researched  to  meet  near-term  goals.  Mid-term  goals 
focus  on  developing  improved  tactical  detection  and  identification  capabilities  for 
both  chemical  and  biological  warfare  agents.  The  focus  of  far-term  efforts  is  on 
multi-agent  sensors  for  CBRN  agent  detection  and  remote/early  warning  CBRN 
detection.  According  to  the  CBRNDP  Annual  Report  to  Congress  (2004),  the  goal  of 
contamination  avoidance  is  “direct  integration  of  CBRN  detectors  as  a  single  system 
into  various  platforms  linked  into  command,  control,  communication,  computer,  and 
intelligence  (C4I)  networks.”  Improvements  in  these  technologies  will  ultimately  lead 
to  better  detection  of  chemical  agents  and  secondary  vapors  contaminating  sites, 
making  it  easier  to  identify  areas  of  the  site  that  are  operating  in  their  optimal  level 
of  MOPP. 

Focus  in  the  area  of  contamination  avoidance  has  several  potential  payoffs. 
CBRN  detection  systems  to  be  developed  will  “provide  the  capability  to  detect, 
identify  in  real  time,  map,  quantify,  and  track  all  known  CBRN  contamination  in 
a  theater  or  operations.”  Acquiring  these  capabilities  will  enable  commanders  of 
a  targeted  site  to  avoid  contamination,  determine  the  need  for  (and  verification 
of)  effective  reconstitution  procedures,  as  well  as  assume  optimal  protection  that 
will  enable  sustainment  of  the  mission  and  continuance  of  fighting  with  minimal 
performance  degradation  and  casualties.  An  important  payoff  of  the  technologies 
being  developed  is  in  the  dual  use  potential.  Occupational  Environment  Health 
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Surveillance  can  use  the  developed  technologies  to  monitor  air  pollution,  noxious 
fumes  inside  enclosed  areas,  and  municipal  water  supplies. 

There  exist  several  major  technical  challenges  for  advancing  contamination 
avoidance  science  and  technology.  These  include  biological  collection,  detection  and 
identification,  improved  agent  discrimination  and  quantification,  sample  processing, 
interferent  (i.e.,  false  positive  and  negative  alarms),  and  ambient  biological  back¬ 
ground  rejection.  Also,  it  is  challenging  to  reduce  size,  weight,  and  the  power  re¬ 
quirement  for  detectors,  as  well  as  power  generation  and  consumption,  development 
of  integrated  biological  and  chemical  detection  systems,  and  fusing  sensor  data  with 
mapping,  imagery,  and  other  data  for  the  near  real-time  display  of  events.  The 
discrimination  capability  of  standoff  detection  for  biological  materials  have  been  im¬ 
proved  in  the  two  years  prior  to  the  publication  of  the  2004  CBRNDP  annual  report 
[5].  Clearly,  the  government  recognizes  the  potential  hostile  use  of  biological  or 
chemical  WMDs  on  the  domestic  United  States  and/or  foreign  interests.  However, 
the  main  focus  of  the  research  has  been  on  avoiding  CBRN  hazards  through  the 
development  of  early  detection  devices  and  a  warning  system.  The  problem  of  this 
thesis  focuses  on  procedures  to  be  employed  in  the  event  of  an  attack  that  could  not 
be  avoided. 

Should  an  attack  occur  on  a  military  installation  or  other  interest,  it  is  im¬ 
portant  to  have  an  effective  response  policy  available.  To  do  this,  the  nature  of  the 
weapons  released  over  the  site  must  be  understood  and  the  behavior  of  the  released 
agent  (e.g.,  biological  or  chemical)  must  be  studied.  Thus,  an  important  aspect 
of  studying  chemical  weapons  attacks  is  understanding  the  spatial  and  temporal 
evolution  of  the  chemical  release.  Several  papers  have  been  written  regarding  emer¬ 
gency  response  subsequent  to  a  biological  or  chemical  attack  [20],  forecast  models 
[17],  source  release  models  [13],  receptor  models,  and  other  topics  concerning  the 
behavior  of  pollution  or  a  released  chemical  [9]. 
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Extensive  research  has  been  conducted  on  anthrax  attacks,  specifically  with 
regard  to  emergency  response  in  populated  areas.  Wein,  Craft,  and  Kaplan  [20] 
compared  various  emergency  responses  to  such  an  attack.  They  develop  a  math¬ 
ematical  model  consisting  of  an  atmospheric  dispersion  model,  an  age-dependent 
dose-response  model,  a  disease  progression  model,  and  a  set  of  spatially  distributed 
two-stage  queueing  systems  consisting  of  antibiotic  distribution  and  hospital  care. 
The  paper  covers  the  response  to  a  biological  weapon  attack,  specifically  anthrax, 
but  not  a  methodology  for  identifying  areas  that  may  become  contaminated.  They 
agree  that  a  Gaussian  plume  model  may  be  “too  simplistic  to  monitor  and  predict 
the  spatiotemporal  anthrax  concentrations”  following  an  attack  but  they  use  the 
“downstream”  portion  with  a  more  sophisticated  atmospheric  model  to  capture  “at¬ 
mospheric  complexities”  that  are  previously  ignored.  Thus,  they  model  the  evolution 
of  anthrax  through  space  and  time  but  do  not  provide  a  methodology  for  searching 
an  area  to  identify  where  people  are  vulnerable  to  inhaling  potentially  deadly  spores. 

Craft,  Wein,  and  Wilkins  [4]  understand  deterrence  is  not  a  reliable  strategy  for 
defense  against  terrorists,  but  the  best  security  against  an  attack  is  in  consequence 
management.  That  is,  how  the  number  of  deaths  can  be  minimized  after  an  attack 
has  occurred.  Again,  the  focus  is  on  treating  individuals  who  suffer  from  anthrax 
inhalation  and  not  on  identifying  a  policy  for  searching  the  area  to  find  where  atmo¬ 
spheric  conditions  have  spread  the  contamination  and  where  time  (or  atmospheric 
condition)  has  reduced  the  risk  of  contamination  to  other  areas. 

Similar  works  have  been  published  for  the  case  of  smallpox.  Kaplan,  Craft, 
and  Wein  [11]  estimate  the  number  of  smallpox  cases  and  deaths  resulting  from  an 
attack  in  a  large  urban  area.  Specifically,  they  show  that  mass  vaccinations  following 
an  attack  result  in  “far  fewer  deaths  and  much  faster  epidemic  eradication”  than  the 
interim  policy  of  isolating  symptomatic  cases  and  performing  traced  vaccinations 
with  mass  vaccinations  as  a  secondary  plan  if  the  first  is  infeasible.  Kaplan,  Craft, 
and  Wein  [12]  also  evaluate  existing  and  alternative  proposals  for  the  emergency 
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response  to  a  smallpox  attack.  Like  the  previously  discussed  papers,  methodology 
for  identifying  contaminated  areas  and  areas  that  are  free  of  contamination  is  not 
considered. 

Developing  an  optimal  sampling  procedure  subsequent  to  a  chemical  agent 
release  over  a  fixed  operational  site,  is  contingent  upon  reliable  spatiotemporal  con¬ 
centration  data.  Research  in  fields  such  as  weather  forecasting,  receptor  models, 
atmospheric  pollution,  vapor  dispersion,  and  other  such  areas  is  necessary  for  de¬ 
vising  ways  to  obtain  this  information.  According  to  Sohev  [17],  a  “demanding” 
application  of  short-term  weather  forecasting  is  the  real-time  emergency  modelling 
of  nuclear  and  chemical  hazardous  events.  The  paper  presents  an  operational  sys¬ 
tem  developed  to  create  a  short-term  forecast  of  potential  risk  areas,  specifically  in 
the  event  of  a  nuclear  power  plant  accident.  The  framework  developed  in  the  work 
displays  the  capability  of  addressing  inverse  problems  of  working  with  an  unknown 
source  of  release.  This  contribution  is  highly  significant  because  the  source  of  a 
chemical  agent  release  is  not  likely  to  be  known. 

Researchers  have  published  substantial  works  on  estimating  the  source  release- 
rate  of  atmospheric  pollution  [14].  Subsequent  to  a  gas  release,  forecasting  the 
concentration  of  the  gas  in  the  atmosphere  depends  on  the  rate  and  the  location  of 
the  release.  Thus,  four  distinct  cases  for  a  single  point  source  of  a  released  gas  have 
been  considered: 

1.  Instantaneous  release  from  a  known  location, 

2.  Instantaneous  release  from  an  unknown  location, 

3.  Extended  release  over  a  period  of  time  from  a  known  location,  and 

4.  Extended  release  over  a  period  of  time  from  an  unknown  location. 

According  to  Kathirgamanathan  et  al.  [14],  the  first  case  is  simple  due  to 
just  a  single  parameter  requiring  estimation.  The  authors  published  a  previous 
work  [13]  for  the  second  case  in  which  they  present  an  inverse  model  to  estimate 
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the  parameters  for  the  model  simultaneously.  Case  three  is  the  focus  of  [14]  where 
the  authors  develop  an  inverse  model  using  methods  from  groundwater  modelling 
literature  and  identify  factors  that  affect  the  accuracy  of  inverse  model  prediction. 
The  authors  state  they  will  address  case  four  in  future  papers. 

A  result  of  the  work  for  case  two  in  [13],  yields  the  advection-diffusion  equation 
for  determining  vapor  concentration  at  specific  points  that  we  refer  to  in  Chapter  3. 
This  equation  is  an  example  of  an  expression  used  for  determining  the  chemical 
concentration  at  specified  points  some  distance  from  the  point  of  release  at  specified 
times.  It  is  important  to  note  that  in  the  event  of  a  chemical  attack,  it  is  unlikely 
that  the  amount  of  chemical  release,  or  the  exact  location  of  the  release,  will  ever  be 
known.  Thus,  even  a  formal  mathematical  model  for  determining  the  concentrations 
at  points  around  the  attack  will  have  only  limited  value  in  discerning  the  actual 
concentration  at  those  areas  unless  a  reading  is  obtained  through  sensors  or  some 
other  instrument. 

Once  methods  for  obtaining  all  data  are  developed,  methods  for  determining 
optimal  sampling  procedures  must  be  developed  to  incorporate  real-time  information 
into  the  methodology.  In  this  thesis,  we  seek  to  identify  an  optimal  sampling  strategy 
for  searching  an  area  subsequent  to  a  chemical  attack.  We  model  the  real  problem  as 
a  network  and  apply  network  analysis  techniques  to  determine  the  optimal  solution. 
Because  network  models  are  so  critical  to  our  approach  to  this  problem,  a  review  of 
rudimentary  concepts  is  provided  next. 

2.3  Network  Modelling  and  Analysis 

Ahuja’s  text  [1]  provides  a  foundation  for  modelling  various  problems  as  net¬ 
work  flow  problems  and  developing  algorithms  to  solve  the  models.  This  text  pro¬ 
vides  the  background  that  is  necessary  for  understanding  the  models  provided  in 
Chapter  3  for  determining  an  optimal  sampling  strategy. 
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A  specific  type  of  network  flow  problem  important  to  this  thesis  is  the  shortest 
path  problem.  According  to  Ahnja  [1],  this  problem  is  “perhaps  the  simplest  of  all 
network  flow  problems.”  However,  that  does  not  diminish  the  importance  of  the 
shortest  path  problem  for  our  purposes.  The  objective  of  the  basic  problem  is  to 
find  a  path  of  minimum  cost  or  length  from  a  source  node  to  a  sink  node  where 
each  has  an  associated  cost  or  length.  We  seek  to  identify  the  shortest  time  path 
while  simultaneously  maximizing  a  reward  (i.e.,  the  number  of  areas  safely  operating 
without  protective  gear). 

The  notation  of  Ahuja  [1]  refers  to  a  directed  network  G  =  (A/”,  A)  where  A f 
is  the  set  of  nodes  and  A  is  the  set  of  arcs  connecting  those  nodes.  There  is  an 
arc  length  or  cost,  ci)3  associated  with  each  (i,  j)  G  A.  The  source  node  s  is  where 
the  flow  (e.g.,  vehicle,  computer  data  packet,  etc.)  originates  and  the  sink  node  t 
is  where  the  flow  terminates.  The  term  A{i )  is  used  to  define  the  arc  adjacency  list 
of  node  i.  That  is  the  list  of  all  arcs  outgoing  from  node  i  into  other  nodes.  The 
“length  of  a  directed  path”  is  defined  to  be  the  sum  of  the  lengths  of  arcs  in  the 
path.  It  is  assumed  the  network  contains  a  directed  path  from  the  source  to  all  other 
nodes  in  the  network,  otherwise,  there  would  exist  isolated  nodes  that  could  never  be 
reached.  Though  there  are  several  types  of  directed  shortest  path  problems,  finding 
the  shortest  path  from  each  node  to  every  other  node  is  the  type  that  is  relevant  to 
this  problem.  It  is  referred  to  as  the  “all-pairs  shortest  path  problem.” 

A  simple  string  model  is  used  to  aid  in  the  understanding  of  a  shortest  path 
problem  between  a  specified  pair  of  nodes  s  and  t.  Ahuja  [1]  points  out  that  this 
model  can  easily  be  extended  to  the  shortest  path  problem  with  multiple  destinations 
and  nonnegative  arc  lengths.  Assuming  the  string  cannot  be  stretched,  use  the  knots 
to  represent  nodes,  where  ct)]  is  the  length  of  string  joining  the  two  knots  i  and  j. 
Taking  hold  of  the  knot  s  in  one  hand  and  the  knot  t  in  the  other  and  pulling  the 
hands  apart  leaves  one  or  more  paths  held  tight.  These  represent  the  shortest  paths 
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from  node  s  to  t  since  they  are  the  lengths  that  are  restricting  the  source  and  sink 
knot  from  being  pulled  farther  apart. 

Several  important  aspects  of  the  shortest  path  problem  can  be  taken  from  this 
string  model  [1].  The  taut  strings  represent  arcs  on  the  shortest  path  and  thus,  the 
distance  between  any  two  successive  nodes  i  and  j  on  that  path  equal  ch]  of  the 
arc  (i,  j)  between  the  nodes.  Considering  any  two  nodes  i  and  j,  that  may  not  be 
successive,  that  are  connected  by  arc  (i,  j)  in  A,  the  shortest  path  distance  from  the 
source  s  plus  ctJ  is  always  as  large  as  the  shortest  path  distance  from  the  source  to 
node  j.  This  composite  distance  may  be  larger  because  the  string  between  knots  i 
and  j  may  not  be  taut  (i.e.  not  in  the  shortest  path). 

According  to  Ahuja  [1],  there  are  two  groups  of  algorithmic  approaches  to 
solving  shortest  path  problems,  both  of  which  are  iterative:  label  setting  and  label 
correcting.  Each  of  these  approaches  assigns  tentative  distance  labels  to  nodes  at 
each  step.  These  distance  labels  are  estimates  of  the  shortest  path  distances.  In 
other  words,  they  are  upper  bounds  on  the  distances.  The  types  of  algorithms  differ 
in  their  methods  of  updating  distance  labels  from  one  step  to  the  next  and  the 
class  of  problems  they  solve.  Label  setting  algorithms  work  by  designating  one  label 
as  permanent  (optimal)  at  each  iteration  and  are  applicable  only  to  shortest  path 
problems  defined  on  acyclic  networks.  The  more  general  label  correcting  algorithms 
consider  all  labels  as  temporary  until  they  are  made  permanent  at  the  final  step. 
They  are  applied  to  all  classes  of  problems  but  are  less  efficient  than  the  label  setting 
algorithm  [1], 

Integer  programming  (IP)  is  used  to  solve  optimization  problems  where  the 
variables  are  discrete  or  are  restricted  to  integer  values.  Specifically,  IP  can  be  used 
to  solve  shortest  path  problems  [21],  The  basic  formulation  has  an  objective  function 
that  seeks  to  maximize  or  minimize  some  profit  or  cost.  For  example,  consider  the 
linear  program  max{cx  :  Ax  <  b,  x  >  0}  where  A  is  an  m  by  n  matrix,  c  an  n- 
dimensional  row  vector,  b  an  m-dimensional  column  vector,  and  x  an  n-dimensional 
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column  vector  of  variables  or  unknowns.  This  becomes  an  integer  program  by  adding 
the  restriction  that  certain  variables  must  be  integers  [21].  There  are  three  different 
forms  of  integer  programming.  A  mixed  integer  program  (MIP)  is  one  where  some 
of  the  variables  must  be  integer  and  others  are  not.  When  all  of  the  variables  are 
integer,  it  is  an  integer  program  (IP)  and  when  all  of  the  variables  are  restricted  to 
values  of  0  or  1,  it  is  referred  to  as  a  binary  integer  program  (BIP)  [21], 

The  Travelling  Salesperson  Problem  (TSP)  is  a  common  binary  integer  pro¬ 
gram  in  operations  research  and  is  very  relevant  to  the  models  presented  in  this 
thesis.  The  original  TSP  is  used  to  model  and  solve  a  problem  where  a  single  sales¬ 
person  must  visit  n  cities  exactly  once  and  return  to  his/her  original  location.  There 
is  a  time  (distance)  qj  associated  with  travelling  from  city  i  to  city  j.  The  objective 
of  the  TSP  is  to  minimize  travel  time  (distance)  by  identifying  the  optimal  order  in 
which  the  salesperson  should  visit  each  city  and  return  to  his/her  starting  point  [21], 
The  TSP  is  directly  related  to  this  thesis  problem  where  the  salesperson  represents 
a  search  crew,  the  cities  represent  the  areas  of  interest,  and  the  time  corresponds  to 
the  time  required  to  travel  between  the  areas.  However,  several  modifications  must 
be  made  to  accommodate  the  restrictions  in  this  thesis  problem  and  these  will  be 
discussed  in  the  following  chapter. 

The  binary  knapsack  problem  is  another  well-known  problem  related  to  this 
thesis.  Essentially,  one  has  a  “knapsack”  and  must  fill  it  with  some  number  of  items. 
Each  item  has  a  cost  (e.g.,  size)  associated  with  it,  and  thus,  using  an  item  takes 
up  some  of  the  resource  (e.g.,  space).  The  objective  is  to  maximize  the  benefit  of 
what  goes  in  the  knapsack  while  adhering  to  the  resource  constraint.  The  knapsack 
problem  relates  to  the  optimal  sampling  problem  because  we  seek  to  maximize  the 
reward  of  the  search  by  choosing  arcs  (with  associated  times)  without  exceeding  the 
available  resource  (time).  However,  the  order  in  which  the  arcs  are  chosen  matters 
to  this  problem  since  they  must  form  a  feasible  path  through  the  system. 
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Several  algorithms  have  been  developed  for  finding  the  least  expected  travel 
time  path  through  a  network.  However,  these  work  only  for  networks  with  deter¬ 
ministic  time-dependent  (or  time- variant)  travel  times  and  for  networks  with  random 
non-time  dependent  (or  time-invariant)  travel  times  [8].  Some  fairly  recent  papers 
have  studied  networks  with  stochastic  and  time- variant  (STV)  travel  times  [16]. 
Miller-Hooks  and  Mahmassani  [16]  state  that  STV  networks  provide  “a  more  appro¬ 
priate  representation”  for  use  in  making  critical  routing  decisions  than  deterministic, 
static  models. 

Miller-Hooks  and  Mahmassani  [16]  present  a  methodology  for  determining 
paths  with  the  least  expected  time  (LET)  through  the  network.  The  resulting  set 
of  paths  are  Pareto-optimal  (or  efficient).  The  results  are  Pareto-optimal  since  no 
other  solution  would  be  optimal  to  all  objectives  in  this  multi-objective  problem. 
The  method  for  determining  LET  paths  for  STV  networks  differs  from  the  method 
for  determining  LET  paths  in  networks  with  time- invariant  travel  times  and  is  more 
difficult.  For  networks  with  random,  time-invariant  travel  times,  each  random  arc 
weight  can  be  set  to  its  expected  value  and  solved  using  deterministic  methods. 

Thomas  and  White  [18]  model  an  anticipatory  route  selection  problem  as  a 
Markov  decision  process  (MDP).  The  objective  of  their  problem  is  to  find  the  “min¬ 
imum  expected  total  cost  route  from  an  origin  to  a  destination  that  anticipates  and 
then  responds  to  service  requests,  if  they  occur,  while  the  vehicle  is  en  route.”  Their 
problem  is  important  to  mobile  communication  technologies  that  enable  communica¬ 
tion  between  dispatchers  and  drivers.  Real-time  information  is  used  to  determine  the 
best  route  for  a  single  pickup  and  delivery  vehicle  when  the  likelihood,  as  a  function 
of  time,  that  a  potential  customer  will  request  a  pickup,  is  known.  The  flexibility 
incorporated  into  the  anticipatory  routing  problem  allows  dispatchers  to  contact  a 
driver,  who  may  be  in  the  process  of  a  pickup  and  delivery,  and  inform  them  of  a 
new  pickup  request.  The  driver  can  then  adjust  his  or  her  route  to  accommodate 
this  new  request.  This  attribute  is  important  to  other  problems  where  the  ability  to 
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receive  real-time  information  allows  the  entity  (e.g.,  truck,  person,  crew)  to  adapt 
its  route  to  accommodate  the  new  information,  thus  optimizing  the  route  [18]. 

The  anticipatory  routing  problem  which  Thomas  and  White  [18]  consider  in¬ 
volves  a  single,  uncapacitated  vehicle  travelling  along  feasible  arcs  from  a  known 
origin  to  a  known  destination  in  a  fixed  amount  of  time.  As  a  vehicle  traverses  an 
arc,  it  accrues  a  travel  cost  dependent  on  the  arc  (e.g.,  time,  money,  resource)  and 
accrues  a  reward  after  visiting  a  node  requiring  a  pickup.  Clearly,  if  no  pickup  is 
requested,  no  reward  is  accrued.  It  is  assumed  each  customer  has  a  known  distri¬ 
bution  giving  the  probability  that  the  customer  will  request  a  pickup  at  time  t  or 
before,  for  all  t.  Thomas  and  White  [18]  formulate  the  model  of  the  road  network  as 
a  MDP.  The  graph  consists  of  a  set  of  nodes  and  edges.  An  arc  exists  if  and  only  if 
there  is  a  road  permitting  travel  between  the  two  endpoints  of  the  arc.  They  define 
a  successor  set  containing  the  adjacent  nodes  for  each  node  which  includes  the  node 
itself  (i.e.,  the  driver  can  wait  at  a  node).  There  is  an  origin  and  destination  node 
in  the  set  of  nodes  but  they  are  not  in  the  set  representing  customers  who  may  or 
may  not  request  service  (i.e.,  requests  cannot  be  made  from  the  origin  or  destination 
node). 

Decision  epochs  for  the  MDP  occur  when  the  vehicle  arrives  at  a  node  and  this 
is  a  finite  horizon  problem  with  the  positive  integer  T  indicating  that  action  selection 
terminates  no  later  than  T  —  1.  A  random  variable  Q  represents  the  total  number 
of  decisions  and  kl(t )  represents  the  state  of  customer  /’s  current  service  request. 
The  successor  set  is  also  the  action  set  whose  cardinality  is  necessarily  finite.  They 
also  define  decision  rules  at  a  time  t  as  a  function  that  selects  an  available  action 
and  a  policy  as  a  sequence  of  those  decision  rules  and  determine  the  structured  cost 
function  and  policy  structure  results.  The  work  is  important  to  this  thesis  because 
it  incorporates  real-time  information  when  making  decisions  by  anticipating  possible 
requests.  The  ability  to  incorporate  real-time  or  near  real-time  information  into 
decisions  allows  for  optimal  policies  being  developed  when  the  problem  is  dynamic. 
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The  main  contribution  of  this  thesis  is  an  optimal  sampling  procedure  for  a 
single  crew  to  assess  the  contamination  level  at  various  areas  throughout  a  fixed 
operational  site  subsequent  to  a  chemical  attack.  The  objective  is  to  identify  the 
maximum  number  of  areas  that  can  safely  operate  without  protective  gear  and  to 
complete  the  search  within  an  allotted  time  limit.  We  consider  four  cases,  increasing 
in  complexity,  and  formulate  each  one  as  a  binary  integer  program  to  develop  a 
methodology  for  optimally  sampling  the  site  in  each  case.  First,  we  consider  static 
and  deterministic  vapor  concentrations  and  then  we  modify  the  problem  to  allow  the 
vapor  concentrations  to  behave  according  to  a  stationary  probability  distribution. 
Next,  we  incorporate  temporal  evolution  into  the  deterministic  model  and  develop 
an  algorithm  for  identifying  the  optimal  sampling  procedure.  Finally,  we  consider 
vapor  concentrations  that  evolve  according  to  a  time- variant  probability  distribution 
and  develop  an  algorithm  that  identifies  the  optimal  sampling  procedure  based  on 
the  probability  of  obtaining  a  reward  at  a  searched  area.  The  final  solution  is  a  time- 
adaptive  policy  directing  the  search  crew  through  the  site  to  the  maximum  number 
of  areas  identified  to  have  vapor  concentrations  below  the  threshold.  Thus,  the 
maximum  number  of  areas  in  which  personnel  can  safely  operate  without  protective 
gear  are  identified. 
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3.  Formal  Mathematical  Model 


In  this  chapter,  we  present  a  mathematical  model  for  optimally  searching  an 
operational  site  subsequent  to  a  chemical  attack.  We  consider  three  distinct  cases 
before  examining  the  fourth  (and  most  realistic)  case,  ultimately  developing  a  policy 
for  searching  the  site,  within  an  allotted  amount  of  time,  to  find  the  maximum 
number  of  areas  in  which  the  level  of  protective  gear  can  be  safely  reduced. 

3.1  Model  Description 

Assume  that  a  fixed  site  (e.g.,  a  military  installation)  is  situated  on  a  Cartesian 
plane  such  that  ordered  pairs  ( x ,  y )  may  be  used  to  identify  and  label  critical  locations 
such  as  the  flight  line,  headquarters,  and  other  operationally  imperative  areas.  In 
Figure  3.1,  such  locations  are  indicated  by  an  asterisk  (*).  Orienting  the  site  on 
the  positive  quadrant  of  the  coordinate  system  with  the  origin  of  the  chemical  agent 
release  in  the  southwest  corner,  each  area  can  be  identified  by  the  ordered  pair 
(x,y),  where  x  and  y  are  the  distances,  in  meters,  from  the  origin  in  their  respective 
directions  to  the  point  where  vapor  concentration  readings  will  be  obtained. 

In  this  research,  we  represent  the  fixed  operational  site  as  a  network  in  which 
the  areas  are  considered  to  be  the  nodes  of  the  network.  Nodes  are  labelled  by  the 


Figure  3.1  Graphical  depiction  of  areas  on  an  installation. 
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integers  1,  2, N,  where  N  is  the  total  number  of  areas.  The  network’s  arcs  connect 
pairs  of  nodes  to  represent  feasible  travel  between  those  nodes.  It  is  important  to 
note,  that  though  the  arcs  are  represented  as  line  segments,  the  actual  path  between 
the  areas  may  not  necessarily  be  straight  (or  direct).  The  arcs  simply  represent 
feasible  travel  paths  from  one  area  to  another. 

Let  G  =  (M,  A),  be  a  directed  graph  where  J\f  is  the  set  of  nodes,  having 
cardinality  N  (i.e. ,  |A/”|  =  N),  and  A  =  {( i,j )  :  i,  j  G  J\f}  is  the  set  of  arcs  connecting 
pairs  of  vertices  i  and  j.  An  arc  (i,j)  is  said  to  be  incident  to  nodes  i  and  j  when  the 
arc  directly  connects  node  i  and  node  j  [1].  A  directed  arc  connecting  nodes  i  and  j 
indicates  feasible  travel  between  nodes  i  and  j  in  the  direction  of  the  arc.  Undirected 
arcs  have  no  arrowheads  on  either  end  and  imply  travel  is  possible  in  both  directions. 
That  is,  one  may  arrive  at  node  j  from  node  i  and  likewise,  one  may  arrive  at  node 
i  from  node  j.  A  leaf  is  a  node  with  at  most  one  arc  incident  thus  allowing  only  one 
route  in  and  the  same  route  out.  Arc  weights  denote  static  and  deterministic  travel 
times  between  two  nodes  and  the  cost  (or  reward)  associated  with  each  node  is  found 
from  the  vapor  concentration  at  that  node.  The  in-degree  of  a  node  is  the  number  of 
arcs  incoming  to  the  node  and  the  out-degree  is  the  number  of  arcs  emanating  from 
the  node.  A  node  j  G  J\f  is  said  to  be  covered  if  it  has  been  visited. 


Figure  3.2  Graphical  depiction  of  network  definitions. 


In  Figure  3.2,  arc  (i,j)  is  incident  to  nodes  i  and  j.  The  deterministic  duration 
tij  is  required  to  travel  from  node  i  to  node  j.  Finally,  at  time  t,  the  vapor  concen¬ 
trations  at  nodes  i  and  j  are  Vi(t)  and  Vj(t),  respectively.  If  the  vapor  concentration 
level  does  not  depend  on  time,  the  notation  is  simply  vt  and  vr 
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The  primary  focus  of  this  thesis  is  the  routing  of  a  search  crew  through  a  fixed 
operational  site  in  order  to  identify  locations  in  which  hazardous  secondary  vapors 
are  present  subsequent  to  a  chemical  weapons  attack.  In  reality,  this  problem  is 
both  dynamic  and  stochastic  in  nature.  That  is,  the  vapor  concentrations  change 
over  time,  either  increasing  in  some  areas  or  decreasing  in  others.  The  problem 
is  stochastic  in  that  the  hazardous  vapor  concentrations  are  random  variables  and 
are  governed  by  some  (time-variant)  probability  distribution  which  is  assumed  to 
be  known.  Due  to  the  stochastic  temporal  and  spatial  evolution  of  the  secondary 
vapors,  it  is  necessary  to  obtain  actual  data  through  chemical  vapor  sensors  and/or 
other  such  instruments.  Technological  capability  to  obtain  such  data  will  contribute 
to  solving  the  problem  by  providing  accurate,  real-time  information.  However,  in 
lieu  of  having  such  data,  we  assume  known  distributions  in  our  model. 

Due  to  the  dynamic  and  stochastic  nature  of  the  vapor  concentration  levels, 
crews  must  intelligently  search  the  site  to  obtain  concentration  measurements  to 
determine  the  correct  level  of  protection  required  for  safe  operations.  Thus,  they 
need  to  reach  as  many  areas  as  possible  within  a  fixed  time  window  while  ignoring 
non-critical  areas  or  areas  not  likely  to  be  contaminated.  In  order  to  approach  this 
problem,  we  first  consider  the  following  three  cases: 

1.  Static  and  deterministic  vapor  concentration  levels; 

2.  Static  and  stochastic  vapor  concentration  levels; 

3.  Dynamic  and  deterministic  vapor  concentration  levels. 

Once  these  three  cases  are  formulated  and  solved,  the  fourth,  and  most  realistic  case, 
of  dynamic  and  stochastic  vapor  concentrations,  can  be  formulated  and  solved  using 
the  techniques  and  insights  obtained  from  the  previous  three  formulations. 

In  the  next  section,  a  number  of  problem  assumptions  are  presented.  Also,  we 
present  the  four  cases  considered  and  the  methodology  for  solving  each  one. 
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3.2  Model  Assumptions  and  Cases 

All  four  cases  employ  the  following  set  of  assumptions: 

1.  The  weapon  released  is  a  chemical  agent  with  known  characteristics  (e.g.,  chem¬ 
ical  makeup  and  properties,  effect  on  human  life,  evaporation  rates,  lifetime, 
etc.). 

2.  Only  one  crew  is  available  to  obtain  instrument  readings  of  vapor  concentration 
levels. 

3.  The  crew  travels  (either  on  foot  or  in  a  vehicle)  at  a  constant  velocity  and 
experiences  zero  delays  due  to  traffic,  etc. 

4.  Travel  times  between  nodes  are  static  and  deterministic. 

5.  Travel  times  between  two  nodes  are  symmetric  (i.e.,  tl3  =  tJti  V(i,  j)  G  -4). 

6.  Only  one  reading  within  each  area  is  required  to  assess  the  correct  level  of 
protective  gear  required. 

7.  The  time  required  to  obtain  an  instrument  reading  is  constant. 

8.  There  is  a  fixed  vapor  concentration  threshold  v*  below  which  protective  gear 
is  not  required  and  above  which  it  is  required. 

9.  Vapor  concentrations  evolve  spatially  (i.e.,  the  level  of  contamination  in  an 
area  depends  on  its  distance  from  the  initial  release). 

The  first  assumption  implies  that  we  focus  on  chemical  weapons,  as  opposed  to 
biological  weapons.  The  next  two  assumptions  regarding  the  search  crew  simplify  the 
problem.  The  methodology  we  develop  is  based  on  a  single  crew  available  for  the  task 
that  travels,  either  in  a  vehicle  or  on  foot,  at  a  constant  velocity.  Again,  to  reduce 
the  number  of  stochastic  elements  in  the  problem,  we  impose  the  assumption  of  a 
fixed  and  deterministic  amount  of  time  required  for  obtaining  instrument  readings. 
This  allows  the  instrument  reading  time  to  be  added  to  the  travel  time.  Also,  we 
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assume  symmetric  travel  times,  thus  implying  travel  from  node  i  to  node  j  requires 
the  same  amount  of  time  as  travelling  from  node  j  to  node  i.  Research  suggests  a 
fixed  vapor  level  exists,  above  which  the  vapor  concentration  is  harmful  to  human 
beings  and  below  which  they  can  operate  safely.  Finally,  we  recognize  the  chemical 
vapor  concentrations  change  as  one  moves  farther  from  the  source  of  release.  This 
assumption  attempts  to  capture  the  true  behavior  of  chemical  vapor  concentrations 
considered  when  developing  the  methodology  for  each  case. 

In  the  following  subsections  we  present  the  three  cases  mentioned  in  section  3.1. 
First,  we  consider  the  case  where  the  vapor  concentration  level  at  each  node  is  static 
and  deterministic.  This  is  a  trivial  case  where  the  vapor  concentrations  can  be  com¬ 
puted  at  the  initial  time  of  the  attack  and  do  not  change  over  time.  Second,  we 
examine  the  case  in  which  vapor  concentrations  are  static  and  stochastic.  This  is 
closer  to  the  real  world  case  because  the  vapor  concentrations  cannot  be  known  with 
certainty  until  an  instrument  reading  is  obtained.  The  probability  distributions  for 
the  vapor  concentrations  are  assumed  to  not  change  over  time.  Third  is  the  case  con¬ 
sidering  dynamic  and  deterministic  vapor  concentrations.  Time  is  incorporated  into 
the  mathematical  formula  and  the  vapor  concentrations  change  (deterministically) 
as  time  progresses.  Finally,  the  fourth  case  considers  dynamic  and  stochastic  vapor 
concentrations.  Representative  of  the  real  world  situation,  the  vapor  concentration 
levels  behave  according  to  a  time-variant  probability  distribution.  There  is  a  clear 
objective  for  each  case  and  the  result  is  a  policy,  in  some  cases  time-adaptive,  that 
directs  a  search  crew  through  the  site  with  the  goal  of  identifying  the  maximum 
number  of  areas  for  which  the  level  of  protective  gear  can  be  safely  reduced. 

3.2.1  Static  and  Deterministic  Vapor  Concentration  Levels 

The  first  case  considers  a  scenario,  subsequent  to  a  chemical  agent  release,  in 
which  the  vapor  concentrations  at  each  area  are  static  and  deterministic  quantities. 
That  is,  they  do  not  change  over  time  and  can  be  determined  via  a  mathematical 
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formula  in  which  all  of  the  parameters  are  known.  An  appropriate  formula,  such  as 
those  vapor  dispersion  models  [14],  can  be  used  to  determine  the  vapor  concentra¬ 
tions,  Vj,  at  each  node.  Figure  3.3  provides  an  example  of  the  network  representation 
for  this  case. 


Figure  3.3  Network  representation  for  a  4-node  site  in  the 
static/deterministic  case. 

This  is  a  multi-objective  case  where  we  seek  to  identify  the  shortest  time  path 
from  the  starting  node  s  to  the  other  nodes  contained  in  the  site  while  maximizing 
the  number  reached.  The  vapor  concentration  levels  are  calculated,  at  a  fixed  time 
to,  from  an  equation  such  as  the  advection-diffusion  Equation  (3.2).  The  following 
parameters  must  be  known  to  employ  Equation  (3.2)  to  compute  Vj,  when  node  j 
has  coordinates  (x,y,  0): 

x,y,z  =  Cartesian  coordinate  system  with  the  rc-axis  oriented  in  the  direction  of 
the  mean  wind,  the  y- axis  in  the  horizontal  cross-wind  direction,  and  the  z-axis 
in  the  upwards  vertical  direction. 

kx,  ky,  kz  =  eddy  cliffusivities  for  the  x,y,  and  z  directions,  respectively  in  m2sec-1 
q  =  the  total  mass  release  in  kg 

h  =  height  above  the  ground  where  the  instantaneous  gas  release  occurs  in  m 
u  =  wind  velocity  in  m/sec 
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The  equation, 
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since  we  assume  z  —  0  (i.e.,  the  x,  y-plane  lays  on  the  ground).  Thus,  from  this  point 
forward,  the  ordered  pair  (x,  y )  is  sufficient  for  denoting  the  coordinates  of  a  node. 


In  the  static  case,  to  denotes  the  elapsed  time  since  the  initial  release  when 
the  vapor  concentration  for  each  node  j  G  J\f  is  calculated  and  it  remains  constant. 
Employing  Equation  (3.2)  yields  the  vapor  concentration  for  coordinates  (x,y)  cor¬ 
responding  to  the  graphical  location  of  each  node  j  G  A/",  with  specified  parameters 
q,  u,  kx,  ky ,  kz,  and  h.  From  these  concentrations,  we  define  a  reward  r3  for  each  node 
j.  For  a  fixed  vapor  concentration  threshold,  denoted  v* ,  the  reward  for  each  node 
is 


{1,  Vj  <  v* 
0,  v3  >  v* 


(3.3) 


The  travel  times  are  used  to  identify  the  optimal  route  the  search  crew  must 
follow.  The  ordered  set  ^  contains  the  nodes  found  in  the  optimal  path  in  the  order 
in  which  those  nodes  should  be  searched.  Following  this  path  allows  the  crew  to 
reach  as  many  areas  as  possible  in  the  minimum  amount  of  time. 

The  mathematical  programming  formulation  for  the  problem  consists  of  a 
multi-objective  function  and  several  constraints.  The  objective  function  (3.4)  seeks 
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to  maximize  the  number  of  areas  sampled  where  a  reward  is  observed  while  also 
minimizing  the  time  required  to  search  the  site.  However,  a  time  constraint  (3.5) 
is  imposed,  where  T  is  the  allotted  time  for  the  search  and  is  the  time  required 
to  travel  from  node  i  to  node  j.  Due  to  (3.5),  the  result  may  be  a  path  that  does 
not  reach  every  node  in  the  network.  The  optimal  path  identified  as  the  solution 
is  denoted,  -0.  In  the  mathematical  formulation,  is  a  binary  variable  assuming 
a  value  of  1  if  the  search  crew  travels  from  node  i  to  node  j  and  0  if  node  j  is  not 
reached  from  node  i.  The  set  S  C  J\f  is  a  set  of  nodes  creating  a  subtour  (i.e.,  a  path 
not  beginning  and  ending  at  the  source  node,  s).  For  networks  whose  nodes  have  a 
minimum  in-degree  of  2,  the  formulation  is 


N  N  N  N 

max  EE  Tjxij(opt ),  min  EE  tijXij 

i=l  j= 1  i=  1  j= 1 


subject  to 


N  N 

i=l  j= 1 

N 

x-^  = 

3= 1 
N 

i=i 

N 

3=1 

xi,j  "F  Xjj  E 


< 


< 


T 

1  for  s  G  J\f 

1  for  j  =  1, ...,  N;j  ±  i 

1  for  i  —  1, ...,  TV;  i  ^  j 
1  for  all  (f,  j)  E  A 


EE*«  <  \S\  -  1  for  5  C  AT,  2  <  \S\  <  TV-  1 
ies  jes 


( 1,  if  (i,j)  E  i\) 
Dj  =  < 

I  0,  otherwise 

Uj  >  0. 


(3.4) 


(3.5) 

(3.6) 

(3.7) 

(3.8) 

(3.9) 

(3.10) 

(3.11) 

(3.12) 
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Clearly,  Equation  (3.6)  requires  that  exactly  one  arc  emanates  from  the  source  node  s 
in  the  final  solution.  This  ensures  the  search  crew  departs  from  its  starting  location. 
The  next  two  equations,  (3.7)  and  (3.8)  restrict  the  number  of  incoming  and  outgoing 
arcs  to  at  most  1.  These  constraints  work  for  the  case  where  the  minimum  in-degree 
of  all  nodes  is  2,  because  the  search  crew  cannot  become  stuck  at  a  leaf  of  the  network 
when  an  outgoing  arc  exists.  Since  there  is  no  benefit  for  returning  to  a  previously 
searched  node  when  the  rewards  are  static,  these  constraints  hold  and  Equation  (3.9) 
is  included  to  prevent  the  crew  from  backtracking  (i.e.,  using  the  same  arc  twice). 
Equation  (3.10)  eliminates  subtours  from  the  solution  as  illustrated  in  Figure  3.4 
and  finally,  (3.12)  restricts  the  times  to  nonnegative  numbers. 


Figure  3.4  Network  representation  of  a  subtour. 


The  mathematical  formulation  bears  similarity  to  the  well-known  Travelling 
Salesperson  Problem  (TSP).  The  objective  function  is  to  minimize  some  time  (dis¬ 
tance)  and  there  is  a  subtour  elimination  constraint.  However,  dne  to  the  knapsack 
problem  constraint  (3.5),  where  time  is  the  “resource,”  it  is  possible  that  not  all 
nodes  will  be  reached,  a  requirement  of  the  TSP.  Also,  the  crew  does  not  have  to 
return  to  the  origin  as  required  in  the  TSP.  Hence,  TSP  solution  procedures  are  not 
directly  appropriate  for  solving  the  problem  in  this  case. 

For  networks  having  a  node  whose  minimum  in-degree  is  1,  it  is  possible  for  the 
crew  to  become  stuck  at  a  leaf  unless  Equations  (3.7),  (3.8),  and  (3.9)  are  removed 
from  the  formulation.  Removing  those  constraints  allows  the  crew  to  backtrack 
along  an  arc  should  it  become  stuck  at  a  leaf  node.  Equations  (3.7)  and  (3.8)  can 
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be  replaced  with 


N 

y:  Xt.j  <  2  for  j  =  1, N;j  ^  i 

i= 1 


^  Xij  <  2  for  j  =  1, N;  i  ±  j, 

3= 1 

respectively.  These  new  constraints  allow  an  arc  to  be  included  at  most  twice, 
allowing  the  crew  to  continue  the  search  from  a  leaf  node. 

We  have  determined  the  best  way  to  solve  this  case  is  by  implicit  enumeration. 
Beginning  with  the  starting  node  s,  each  adjacent  arc  is  examined.  If  the  current 
time  t  plus  the  next  travel  time  is  greater  than  the  allotted  time  T,  the  branch 
terminates  at  node  i.  Due  to  this  time  constraint,  an  upper  bound  can  be  found  on 
the  number  of  possible  arcs  in  the  optimal  solution.  This  upper  bound  is  given  by 


u.b. 


T 


(3.13) 


Determining  this  upper  bound  reduces  the  size  of  the  problem  by  providing  a  limit 
to  the  number  of  arcs  in  the  optimal  path.  For  example,  if  T  =  60  minutes  is  allotted 
for  the  search,  and  the  minimum  travel  time  for  all  arcs  was  5  minutes,  the  upper 
bound  on  the  number  of  arcs  in  the  optimal  solution  is  12.  Hence,  if  there  were 
more  than  12  nodes  in  the  network,  not  all  of  them  could  be  reached  in  the  optimal 
solution.  Also,  the  size  of  the  problem  is  further  reduced  by  restricting  backtracking 
unless  necessary  and  by  not  allowing  the  same  node  to  be  searched  more  than  once. 

We  define  a  branch  as  an  ordered  vector  of  nodes  representing  a  feasible  path 
with  an  associated  time  length,  denoted  Ty.  Each  branch  6*.  for  k  —  1,  where 
B  is  the  number  of  branches,  is  examined  and  if  the  cardinality  of  the  branch  equals 
the  number  of  nodes  in  the  network,  the  path  of  minimum  time  is  chosen.  That 
is,  for  all  branches  such  that  \bk\  =  N,  the  optimal  path  x[)  is  the  branch  b*k  with 
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minimal  total  search  time 


Tk  =  min  {rk}  (3.14) 

If  the  maximum  cardinality  of  all  branches 

b  =  max{|6fc|}  (3.15) 

is  less  than  number  of  nodes,  the  branch  of  maximum  cardinality  (and  shortest  time) 
is  optimal.  That  is,  the  optimal  path  0  =  b*k  is  that  path  for  which 

tI  =  min  {by  :  |6fc|  =  b}.  (3.16) 

k 

Upon  completion  of  the  path  specihed  by  the  optimal  path  0  with  total  time 
70,  instrument  readings  are  obtained  for  all  nodes  in  the  path  (i.e.,  Vj  G  0).  The 
term  r3 ,  represents  the  indicator  function,  I{vj  <  v*},  and  is  1  if  Vj  <  v*  and  0 
otherwise.  That  is,  if  the  vapor  concentration  at  node  j  is,  in  fact,  less  than  the 
predetermined  threshold,  then  the  indicator  function  will  assume  the  value  unity. 
For  all  nodes  not  in  the  path,  j  0  y the  vapor  concentration  v3  obtained  from 
Equation  (3.2)  is  used  to  determine  the  reward  value  r3  for  those  nodes.  The  total 
reward  r*  is  the  sum  of  all  reward  values,  denoted  by 

r’  =  J2r>-  <3'17) 

jeM 

Next  we  consider  the  case  in  which  the  vapor  concentrations  are  assumed  to 
behave  according  to  a  known  probability  distribution  that  is  time-invariant. 

3.2.2  Static  and  Stochastic  Vapor  Concentration  Levels 

The  second  case  considers  a  scenario  in  which  the  vapor  concentration  at  each 
node  is  distributed  according  to  some  known  (and  time-invariant)  probability  dis- 
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tribution  with  a  finite  first  moment.  That  is,  the  vapor  concentration  level  at  node 
j  is  not  known,  but  rather  the  expected  value  of  the  vapor  concentration  at  node  j 
can  be  computed  or  estimated.  The  random  variable  V)  denotes  the  random  vapor 
concentration  at  node  j. 

It  has  been  shown  by  Hall  [8]  that  the  least  expected  travel  time  path  through 
time-invariant  stochastic  networks  can  be  solved  using  the  same  methods  as  those 
for  time-invariant  deterministic  networks.  Thus,  in  this  static  case  the  problem 
can  be  solved  using  implicit  enumeration.  Figure  3.5  is  an  example  of  the  network 
representation  of  this  case. 


E[V,]  E[V2] 


E[V3]  E[V4] 

Figure  3.5  Network  representation  of  a  4-node  site  in  the 
static/stochastic  case. 

We  seek  to  maximize  the  reward  while  minimizing  the  time  required  to  travel 
through  the  network.  This  allows  the  crew  to  reach  as  many  nodes  as  possible 
in  the  allotted  time.  The  difference  from  the  previous  case  is,  before  the  search 
begins,  vapor  concentrations  at  each  node  j  cannot  be  determined  with  an  equation 
such  as  Equation  (3.1),  but  the  expected  value  can  be  computed  through  a  known 
(stationary)  probability  distribution.  Thus,  using  the  expected  vapor  concentration, 
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the  reward  value  for  each  node  is  determined  by 


Jl,  E[V,\<v- 
\o,  E\Vj]  >  v‘ 


(3.18) 


Here,  expected  vapor  concentrations  are  considered  rather  than  deterministic  values 
and  actual  vapor  concentrations  cannot  be  known  until  instrument  readings  are 
obtained. 

Unlike  the  case  in  section  3.2.1,  the  initial  rj  values  are  based  on  the  expected 
vapor  concentrations.  Therefore,  a  realization  of  the  vapor  concentration  as  a  result 
of  the  crew’s  search  may  reveal  more  areas  below  the  threshold,  or  more  areas  above 
the  threshold,  than  initially  expected.  Essentially,  the  true  Vj  values  cannot  be 
assigned  to  an  area  until  the  crew  obtains  the  instrument  reading  for  that  area.  Thus, 
the  importance  of  searching  as  many  areas  in  the  allotted  time  is  more  notable  for 
this  case  where  the  actual  total  reward  can  only  be  known  following  an  instrument 
reading  of  every  site.  At  time  T,  the  end  of  the  time  allotted  for  the  search,  the 
decision  makers  will  know  the  total  reward  r*  such  that 
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In  the  next  section  we  consider  the  case  where  the  vapor  concentrations  are 
deterministic,  but  dynamic.  We  present  an  algorithm  for  determining  the  optimal 
policy  that  will  maximize  the  reward  while  completing  within  the  allotted  time. 


3.2.3  Dynamic  and  Deterministic  Vapor  Concentration  Levels 

The  third  case  is  distinct  from  the  Erst  two  cases  in  that  we  now  assume  the 
vapor  concentrations  depend  on  time.  The  vapor  concentration  for  each  node  can 
presumably  be  computed  from  known  physical  relationships  but  the  values  depend 
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on  the  time  of  the  calculation.  This  is  more  representative  of  the  real  world  case 
since  chemical  vapor  concentrations  do  change  over  time  (and  space).  Figure  3.6 
provides  an  example  of  the  network  representation  for  this  case. 


Vi(t)  V2(t) 


Figure  3.6  Network  representation  of  a  4-node  site  in  the  dy¬ 
namic/deterministic  case. 


As  before,  let  Vj(t)  denote  the  vapor  concentration  at  node  j  at  time  t.  The 
advection-diffusion  Equation  (3.2)  is  an  example  of  how  one  might  calculate  v3{t) 
where  node  j  is  located  at  coordinate  (x,  y ).  In  this  case,  we  allow  t  to  vary,  allowing 
the  vapor  concentration  for  each  node  j  G  J\f  to  be  updated  at  various  times.  Recall 
the  equation 


Vj(t) 


Q 

— ^ - exp 

87T2  (kxkykzy/2t3/2 


(x  —  ut )2 


Akyt 


(3.19) 


where  t  is  updated  to  represent  the  time  the  crew  will  arrive  at  the  next  node  j  (i.e. , 
t  —  t  +  Uj  when  the  current  location  is  at  node  i).  Using  this  equation,  the  search 
crew  can  calculate  the  vapor  concentrations  for  all  areas  adjacent  to  their  current 
location  at  the  time  of  their  arrival  to  that  subsequent  area. 

Again,  the  decision  maker’s  desire  is  to  return  the  site  to  full  operational  capa¬ 
bilities  as  soon  as  possible.  Thus,  they  need  to  know  where  and  when  the  protective 
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gear  requirement  can  be  withdrawn  and  where  and  when  it  needs  to  be  implemented. 
Therefore,  we  develop  an  algorithm  that  seeks  to  maximize  the  number  of  areas 
searched  in  which  the  protective  gear  is  no  longer  required  (i.e. ,  where  the  reward 
value  is  1). 

The  algorithm  begins  with  an  initialization  step  that  commences  at  to,  the  time 
duration  from  when  the  chemical  agent  was  released  over  the  site  to  the  completion 
of  the  first  instrument  reading.  The  crew  uses  an  appropriate  instrument  for  reading 
vapor  concentrations  to  determine  the  concentration  for  the  origin  node  s.  Thus, 
the  reward  for  node  s  is 


rs(t o) 


j  1,  Vj(t0)  <  V* 

[o,  Vj(t0)>v* 


(3.20) 


The  algorithm  then  considers  all  areas  adjacent  to  area  s.  Vapor  concentrations  are 
calculated  for  those  areas  using  Equation  (3.19)  and  the  time  it  would  be  when  the 
crew  arrives  there  (i.e.  the  current  time,  to,  plus  the  travel  time  from  s  to  j,  to  +  tsj). 

If  an  area  has  a  future  vapor  concentration  level  below  the  threshold,  it  is  a 
candidate  for  reduction.  Each  such  area  is  then  considered  individually  and  the  one 
with  the  minimum  difference  from  the  threshold  is  chosen  (i.e.,  the  node  such  that 
Vj(t)  is  closest  to  v*).  We  choose  the  area  with  the  minimum  difference  below  the 
threshold  because  it  is  the  most  critical  area.  That  is,  it  is  one  that  most  demands  an 
instrument  reading  since  it  is  closest  to  being  considered  contaminated  to  a  degree 
that  is  unsafe  for  personnel  working  in  the  area  without  protective  gear.  If  there  are 
no  areas  below  the  threshold,  the  node  within  the  shortest  amount  of  time  (distance) 
is  chosen.  This  is  because  the  crew  must  leave  its  current  location  and  if  no  adjacent 
node  has  a  positive  reward,  choosing  the  closest  node  will  advance  the  crew  in  a 
minimum  amount  of  time  so  more  nodes  can  be  considered. 
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This  method  is  repeated  for  all  adjacent  areas  and  the  time  t  is  updated  ap¬ 
propriately  until  the  allotted  search  time  has  expired.  Previously  searched  nodes 
can  be  revisited  if  the  vapor  concentration  at  a  later  time  is  found  to  be  below  the 
threshold.  As  measurements  are  taken  in  each  area,  the  value  of  r,-  is  assigned.  When 
the  search  time  expires,  decision  makers  will  know  the  exact  number  and  location 
of  areas  operating  safely  without  protective  gear  and  those  still  requiring  the  gear 
for  safe  operations.  The  formal  algorithm  for  solving  the  mathematical  program  is 
provided  in  what  follows.  The  set  A/"  denotes  the  set  of  N  nodes  in  the  area  to  be 
searched  and  J\ft  is  the  set  of  nodes  adjacent  to  node  i  G  J\f.  The  set  1Z  is  used  to 
record  nodes  in  the  path  with  a  reward  value  of  1  and  the  set  ^  is  the  ordered  set  of 
nodes  in  the  optimal  path. 

Initialization: 


{1,2,  ...,7V}; 

M  =  O'  :  (i,j)  e  A}; 
ft  =  0; 

^  =  0; 
t  <  t0, 
i  =  1; 

Calculate  current  vapor  level  at  node  z,  Vi(t) 
If  Vi(t)  <  v* 

n(t)  <- 1; 

n  =  {i}-, 

^  <-  ^  U  {i}; 

Else 

Ti(t)  0; 
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^  <-  ^  U  {i}; 

End 
Step  1 

Calculate  Vj(t  +  ttj)  Vj  G  J\ft 
If  Vj(t  +  tij)  <  v* 
rj(t  +  Uj)  *-  l; 

Else 

rj(t  +  Uj)  *-  0; 

End 

Step  2 

For  each  j  G  J\f  such  that  r3  =  1 

Choose  j  such  that  Vj(t  +  tl3)  =  arg rnirg^v; {u*  —  v3(t  +  Cj)} 

n  <-  7eu{j}; 

^  <-  ^  U  {j}; 

End 

if  r.,  (f  +  tij)  =  0  V  j  G  A4: 

Choose  j  such  that  ttJ  =  mug  {tj  j } Vj  G  A/* 

^  <-  ^  U  {j}; 

End 

t  <  t  T  tij , 

Step  3 

If  t  >  T 

STOP 

Else 
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Return  to  Step  1 


End 

This  algorithm  has  been  coded  in  the  standard  mathematical  computing  en¬ 
vironment  MATLAB®  to  determine  a  policy  for  searching  the  site  when  the  vapor 
concentration  levels  are  assumed  to  be  dynamic  and  deterministic.  The  order  in 
which  the  nodes  are  added  to  the  set  ip  is  the  order  in  which  the  site  should  be 
searched.  Moreover,  the  cardinality  of  the  set  1Z  is  the  number  of  areas  where  the 
protective  gear  is  no  longer  required  for  the  safety  of  the  personnel  working  in  those 
areas.  Ultimately,  the  algorithm  results  in  a  time-adaptive  policy  that  directs  the 
crew  to  the  maximum  number  of  areas  where  it  is  critical  to  search  while  reaching 
as  many  as  possible  within  the  allotted  time. 

The  next  subsection  considers  the  final  and  most  realistic  case  in  which  the 
vapor  concentration  levels  are  both  dynamic  and  stochastic.  This  more  complicated 
case  is  approached  in  a  manner  similar  to  the  one  used  in  the  previous  case. 

3.2.4  Dynamic  and  Stochastic  Vapor  Concentrations 

Finally,  the  fourth  case  is  most  representative  of  the  real  world  problem  where 
the  vapor  concentrations  at  each  node  are  denoted  by  time-dependent  random  vari¬ 
ables.  The  solution  to  this  final  problem  will  also  be  a  time-adaptive  policy. 

If  a  chemical  agent  is  released  over  an  area,  the  hazardous  secondary  vapors 
will  continue  to  contaminate  other  areas  and  pose  harm  to  the  personnel  working 
in  the  those  areas.  This  case  partially  captures  the  true  nature  of  those  secondary 
vapors  by  not  only  incorporating  time  into  the  formulation,  but  also  assuming  a 
time-variant  probability  distribution  to  determine  the  probability  that  the  vapor 
concentration  at  each  area  of  interest  is  below  the  threshold  level. 
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Figure  3.7  illustrates  a  possible  network  representation  for  this  case.  The  label 
on  each  node,  Vj(t)  denotes  the  random  vapor  concentration  level  at  node  j  at  time 
t  and,  as  before,  each  arc  weight  represents  the  travel  time  from  node  i  to  node  j. 


v,(t)  V2(t) 


V3(t)  v4(t) 


Figure  3.7  Network  representation  of  a  4-node  site  in  the  dy¬ 
namic/stochastic  case. 

We  assume  each  node  has  an  associated  probability  distribution  that  governs 
the  behavior  of  the  vapor  concentration,  over  time,  at  that  node.  As  time  progresses, 
the  parameters  of  the  distribution  change.  In  this  work,  we  assume  the  concentration 
at  each  node  is  a  random  variable  with  an  exponential  distribution.  Each  node  j 
will  have  an  associated  rate  parameter  at  time  t,  namely,  t-i3{i).  The  cumulative 
distribution  function  for  Vj(t)  is 


Tljit)  =  P(Vj(t)  <  v) 


{1  —  exp ■  v)  v  >  0 
0  v  <  0 


(3.21) 


Suppose  the  rate  parameter  for  node  1  at  t  —  420  sec  is  //i(420)  =  144.30.  The 
probability  that  node  1  is  less  than  the  threshold  v*  =  0.0006 mg/m3  at  time  420 
seconds  is  0.0829.  This  probability  is  calculated  by 


7^(420)  =  1  -  exp (—144.30  •  0.0006)  =  0.0829.  (3.22) 
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Continuing  in  this  fashion,  the  probability  7r ,(£)  that  each  node  is  less  than  the 
threshold  can  be  determined  for  each  time  t.  These  probability  values  are  used  in 
the  algorithm  presented  in  this  section. 

The  algorithm  we  developed  for  this  case  has  a  structure  similar  to  that  of  the 
case  assuming  deterministic  and  dynamic  vapor  concentrations.  It  begins  with  an 
initialization  step  at  time  to  (the  duration  of  time  since  the  chemical  agent  release 
and  the  completion  of  the  first  instrument  reading),  and  with  the  search  crew  located 
at  node  s  (the  starting  point  for  the  search).  The  crew  obtains  an  instrument  reading 
for  its  starting  node  and  each  node  it  searches.  If  Vj(t)  is  below  the  threshold,  the 
node  receives  a  reward  value  at  time  rj(t)  according  to 


fj{t) 


1 1,  Vj{t)  <  v* 
I  0,  Vj(t)  >  v* 


If  r.j  =  1,  the  set  of  reduced  nodes  1Z  is  updated  by 


n  =  nu{j}, 


(3.23) 


(3.24) 


indicating  the  vapor  concentration  at  node  j  has  dropped  below  the  threshold  for 
time  t. 

For  each  adjacent  node,  the  cumulative  distribution  function  (c.d.f.)  of  the 
appropriate  probability  distribution  is  used  to  determine  the  probability  that  the 
next  node  to  be  searched  will  have  a  vapor  concentration  level  below  the  threshold 
upon  arrival.  The  crew  then  chooses  the  area  with  maximum  probability  of  having 
a  reward.  This  will  direct  them  to  the  area  where  they  are  most  likely  to  obtain  a 
vapor  concentration  reading  below  the  threshold  at  the  next  time  step.  Once  the 
crew  arrives  at  the  area,  it  must  obtain  an  instrument  reading  to  determine  the  actual 
vapor  concentration.  In  the  algorithm,  this  is  accomplished  by  drawing  a  value  from 
the  probability  distribution  as  the  realization  Vj{t)  of  the  random  variable  Vj(t).  If 
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the  vapor  concentration  is  below  the  threshold,  the  reward  value  for  that  node,  at 
time  t,  assumes  a  value  of  1  and  is  added  to  the  set  7 Z. 

There  is  a  single  stopping  criterion  terminating  the  algorithm.  Once  the  time 
allotted  for  the  search  expires  the  search  crew  must  cease  searching  the  site  and 
end  at  its  current  location.  Upon  termination,  the  set  7 Z  contains,  in  order,  the 
nodes  searched  where  the  vapor  concentration  was  decreased.  The  set  0  contains 
the  nodes,  in  order,  to  be  searched  within  the  allotted  time  to  maximize  the  reward. 

For  nodes  that  remain  uncovered,  the  expected  value  of  the  vapor  concentration 
is  used  to  determine  the  reward  value.  For  uncovered  nodes,  the  reward  is 


1,  E[Vj(t)]  <  v • 
jo,  E[V,(t)]>v* 


(3.25) 


Thus,  the  final  reward  is 

r*(r  *)  =  5>,W,  (3.26) 

jeM 

such  that  t  <T  is  the  time  of  the  last  instrument  reading  and  r*  is  the  time  length 
of  the  solution.  This  value,  r*(r*)  is  the  total  number  of  nodes  that,  according  to 
actual  instrument  readings  for  j  G  0  and  the  rewards  for  j  ^  0  (determined  from 
the  expected  vapor  concentrations),  no  longer  require  protective  gear. 

Upon  completion  of  the  algorithm,  decision  makers  will  know  which  areas  of  the 
site  are  operating  safely  without  protective  gear,  those  that  still  require  the  gear  for 
safe  operations,  those  which  are  expected  to  not  require  the  gear,  and  those  expected 
to  still  require  the  gear.  Below  is  the  formal  algorithm  for  solving  the  mathematical 
program.  Recall,  the  set  J\f  denotes  the  set  of  N  nodes  in  the  area  to  be  searched 
and  TV)  is  the  set  of  nodes  adjacent  to  node  i  G  TV”.  The  set  1Z  is  used  to  record  nodes 
in  the  path  with  a  reward  value  of  1  and  the  set  0  is  the  ordered  set  of  nodes  in  the 
optimal  path. 


3-21 


Initialization: 


AT  =  {1,2, 

M  =  O'  :  (i,j)  e  .A}; 

72  =  0; 

^  =  0; 
t  <—  t0; 

*  =  i; 

Obtain  current  vapor  level  Uj(t) 

If  Vi(t)  <  v* 

n(t )  <-  l; 

72.  <—  7?.  U  {i}; 

^  <-  ^  U  {i}; 

Else 

Ti(t)  <-  0; 

^  <-  ^  U  {i}; 

End 
Step  1 

Calculate  7r,-(t  +  Cj)  =  P{Vj(t  +  th])  <  u*}  Mj  e  TV; 

Step  2 

Choose  j  such  that  vr^t)  =  max^g^  P{Vj(t  +  titj )  <  u*} 
Obtain  instrument  reading  at  this  node. 

If  Vj(t  +  tjj)  <  u* 
rjit)  <-  1; 
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n  <-  nu{j }; 

^  <-  ^  U  {j}; 
t  ^  +  tij  5 

Else 

rj(t)  <-  0; 

^  <-  ^  U  {j}; 
t  <  t  T  tj  j , 

End 

Step  3 

If  t  >  T 

STOP 

Else 

*  <-  i; 

Return  to  Step  1 

End 

This  algorithm  differs  from  the  previous  algorithm  due  to  the  stochastic  nature 
of  the  vapor  concentration  levels.  It  cannot  be  known  with  certainty  whether  or  not 
an  area  has  a  vapor  concentration  below  the  threshold  until  it  is  searched  by  the 
crew  and  an  instrument  measurement  obtained.  Thus,  the  algorithm  chooses  which 
nodes  to  search  based  solely  on  the  probability  the  vapor  concentration  level  will  be 
below  the  threshold  upon  the  crew’s  arrival. 

Since  there  exists  a  non- zero  probability  that  a  searched  area  will  have  a  vapor 
concentration  above  the  threshold,  the  algorithm  allows  the  crew  to  search  that  area 
again  by  not  removing  it  from  the  set  of  candidate  nodes.  The  resulting  policy  can 
include  searching  a  node  more  than  once  since  the  probability  distributions  change 
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over  time,  and  the  nodes  are  searched  based  on  probabilities  which  may  not  be  very 

high. 

The  next  chapter  provides  several  numerical  examples  to  illustrate  each  of 
these  four  cases.  We  setup  each  case  and  employ  the  mathematical  formulations  and 
algorithms  developed  in  this  chapter  to  find  an  optimal  sampling  procedure  for  the 
site. 
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4.  Numerical  Results 


In  this  chapter,  the  main  results  of  Chapter  3  will  be  illustrated  through 
four  example  problems.  The  algorithms  developed  in  the  previous  chapter  will  be 
used  to  identify  a  policy  that  directs  the  search  crew  through  the  site  in  order  to 
identify  the  maximum  number  of  areas  at  which  protective  gear  may  be  removed.  We 
will  compare  the  results  for  the  static/deterministic  with  the  dynamic/deterministic 
vapor  concentration  levels,  as  well  as  the  results  for  the  static/stochastic  with  the 
dynamic/stochastic  vapor  concentration  levels  to  discuss  the  impact  of  excluding  the 
dynamic  nature  of  the  secondary  vapor  concentrations  when  searching  the  site. 

4-1  Generating  Problem  Instances 

We  randomly  generate  ten  ordered  pairs,  {(ay,  jy)  :  i  =  1,2,...,  10},  to  represent 
a  site  with  ten  critical  areas.  Plotting  these  coordinates  depicts  the  location  of  the 
critical  areas  with  respect  to  each  other  and  the  initial  chemical  agent  release.  From 
this  plot,  a  network  graph  is  drawn  using  nodes  to  represent  the  critical  areas  and 
arcs  to  represent  feasible  travel  between  pairs  of  areas. 

We  arbitrarily  consider  a  rectangular  area  of  the  site  100  meters  by  150  me¬ 
ters.  A  standard  spreadsheet  application  was  used  to  randomly  generate  the  ten  x- 
coordinates  and  ten  ^-coordinates  recorded  in  Table  4.1.  The  points  are  distributed 
uniformly  over  the  intervals  (5, 100)  and  (5, 150),  respectively,  and  the  pairs  are  given 
in  ascending  order  by  the  x  value  to  assign  node  numbers  to  each  area. 

This  example  assumes  all  critical  areas  of  the  site  lie  to  the  northeast  of  the 
chemical  agent  release.  If  the  release  occurred  elsewhere,  the  origin  may  be  translated 
to  the  appropriate  location  without  affecting  the  solution  procedure.  The  areas  are 
plotted  in  Figure  4.1. 
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Table  4.1  x,  ^-coordinates  for  example  site. 


Node 

x  (m) 

y  M 

1 

6.276 

23.015 

2 

14.466 

56.527 

3 

33.297 

96.017 

4 

42.992 

5.234 

5 

60.240 

53.504 

6 

71.900 

140.526 

7 

86.365 

36.667 

8 

94.149 

86.361 

9 

98.829 

7.664 

10 

99.742 

100.500 

Creating  a  node  for  each  point  on  the  plot  in  Figure  4.1,  and  connecting  pairs 
of  nodes  with  arcs,  generates  the  network  representation  for  this  example.  Also 
assuming  zero  detours,  the  Euclidean  distance  from  node  i  to  node  j,  denoted  dt.j, 
can  be  computed  as 

dij  =  \] ( xi  ~  xjY  +  (Vi  ~  Vj)2-  (4-1) 

Assuming  zero  delays  or  obstacles,  and  a  constant  velocity  of  26.8224  m/min  (1 
rni/h)  we  calculate  the  travel  time  for  a  search  crew  from  one  node  to  another.  Five 
minutes  are  added  to  the  time  on  each  arc  (i,j)  G  A  to  account  for  the  instrument 
reading  at  node  j.  Recall,  the  travel  times  are  symmetric  so  that  thj  =  tJ  t  for  all 

(*,  j)  e  A. 

Figure  4.2  depicts  the  network  graph  for  the  notional  site.  The  following  sec¬ 
tions  use  this  site  as  a  template  for  each  example.  The  appropriate  methodology  is 
applied  in  each  section  to  obtain  a  path  directing  the  search  crew  through  the  site 
to  identify  the  maximum  number  of  areas  at  which  the  protective  posture  level  can 
be  reduced. 
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Figure  4.1  Plot  of  x,  ^/-coordinates  for  the  notional  site. 

4-2  Example  1:  Static  and  Deterministic  Vapor  Concentrations 

Recall  for  this  case  that  vapor  concentrations  are  time-invariant  and  can  be  cal¬ 
culated  using  a  deterministic  formula  such  as  the  advection-diffusion  Equation  (3.1) 
(see  Chapter  3).  We  seek  a  policy  ^  to  identify  the  areas  of  the  site  where  the  vapor 
concentration  has  decreased  below  the  threshold  since  the  initial  chemical  agent  re¬ 
lease.  A  decrease  below  the  threshold  value  v*  allows  personnel  in  the  area  to  carry 
out  the  mission  safely  without  the  cumbersome  protective  gear.  The  search  crew  will 
use  this  policy  to  search  the  site  and  obtain  instrument  readings  for  those  areas. 

We  begin  by  calculating  the  vapor  concentration  level  for  each  node  at  initial 
time  to  (i.e.,  time  since  the  chemical  agent  release).  Parameters  used  in  Equa¬ 
tion  (3.1)  are  tabulated  in  Table  4.3. 

The  calculated  vapor  concentration  and  reward  value  at  each  node  is  recorded 
in  Table  4.4.  The  reward  value  is  calculated  as  follows  for  threshold  v*  =  0.0006  mg/m3: 


j  !,  vj  <  v*, 
I  0,  Vj  >  v* 


(4.2) 
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Table  4.2  Distance  and  travel  time  between  nodes. 


Arc 

Lengthen) 

Travel  Time  (min) 

(1.2) 

34.50 

6.29 

(1.4) 

40.79 

6.52 

(1.5) 

61.98 

7.31 

(2.3) 

43.75 

6.63 

(2,5) 

45.87 

6.71 

(3,4) 

91.30 

8.40 

(3,5) 

50.33 

6.88 

(3,6) 

58.92 

7.20 

(4,7) 

53.56 

7.00 

(4,10) 

110.89 

9.13 

(5,7) 

31.08 

6.16 

(5,8) 

47.22 

6.76 

(6,8) 

58.56 

7.18 

(6,9) 

135.56 

10.05 

(7,10) 

65.22 

7.43 

(8,9) 

78.84 

7.94 

(8,10) 

15.20 

5.57 

The  travel  times,  in  minutes,  are  calculated  using  the  relationship 


til  =  <^  +  300, 
v 


(4.3) 


where  v  denotes  the  constant  velocity  of  the  search  crew  (in  m/sec)  and  five  minutes 
(300  sec)  required  for  the  instrument  reading  is  added  to  the  travel  time.  These 
times  are  entered  into  the  matrix  T  =  [Uj],  where  a  dash  (-)  for  th]  signifies  travel 
is  infeasible  between  nodes  i  and  j.  The  times  are  symmetric  and  remain  constant 
throughout  the  duration  of  the  search. 
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Figure  4.2  Network  representation  for  the  10-node  network. 
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The  optimal  path  for  this  example  problem  was  obtained  using  code  written  in 
the  high-performance  language  MATLAB®.  The  code  implicitly  enumerates  every 
feasible  path  that  the  search  crew  can  complete  in  the  allotted  time  of  T  =  90 
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Table  4.3  Parameter  values  for  static/deterministic  example. 


Parameter 

Description 

Value 

to 

time  since  release 

420  sec 

q 

amount  of  chemical  agent  released 

1000  kg 

eddy  diffusivity 

12  m2/sec 

ky 

eddy  diffusivity 

12  m2/sec 

kz 

eddy  diffusivity 

0.2113  m2/sec 

u 

wind  velocity 

0.25  m/sec 

h 

height  of  release 

50  m 

V 

vapor  concentration  threshold 

0.0006  mg/m3 

Table  4.4 


Vapor  concentration  and  reward  value  for  each  node  in  the 
static/deterministic  example  (v*  =  6.0  x  10-4). 


Node 

Vj  (x  10"3) 

ri 

1 

0.496 

1 

2 

0.470 

1 

3 

0.405 

1 

4 

0.682 

0 

5 

0.649 

0 

6 

0.294 

1 

7 

0.760 

0 

8 

0.567 

1 

9 

0.822 

0 

10 

0.500 

1 

minutes  to  identify  the  policy  ijj  that  maximizes 
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ieAf  j&N 

The  solution  to  this  problem,  where  we  arbitrarily  chose  node  s  =  5  as  the  starting 
location  is  illustrated  in  Figure  4.3.  This  solution  yields  the  integer  reward  of  6 
which  is  the  total  number  of  areas  whose  vapor  concentration  is  below  the  threshold 
level  of  v*  =  0.0006  mg/m3.  The  path 


^  =  {5,8,10,7,4,1,2,3,6,9} 
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directs  the  search  crew  through  the  network  and  is  able  to  reach  each  node  exactly 
once  in  a  total  time  of  62.85  minutes.  There  is  a  positive  reward  of  1  received  from 
searching  each  of  the  nodes  in  the  set 

77  =  {1,2,3,6,8,10}. 

No  other  path  exists  that  reaches  all  ten  nodes  in  less  time.  More  importantly, 
the  path  -0  directs  the  crew  through  the  site  and  reaches  each  of  the  critical  areas 
where  the  level  of  MOPP  can  be  reduced.  The  bold  nodes  in  Figure  4.3  represent 
those  areas  of  the  site  where  protective  gear  is  no  longer  required  and  the  numbers 
on  the  arcs  represent  their  order  in  the  search.  We  provide  an  illustration  for  the 
static  and  stochastic  case  in  the  next  section. 


Figure  4.3  Optimal  sampling  path  when  concentrations  are  static 
and  deterministic. 
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4-3  Example  2:  Static  and  Stochastic  Vapor  Concentrations 


For  the  case  of  static  and  stochastic  vapor  concentrations,  the  vapor  concentra¬ 
tion  level  at  each  node  cannot  be  determined  via  a  deterministic  formula.  Instead,  we 
assume  the  concentration  level  is  a  time-invariant  and  continuous  random  variable. 

We  determine  the  expected  value  of  the  vapor  concentration  for  each  node  at 
time  to-  Thus,  time-invariant  probability  distributions  and  appropriate  parameters 
are  assigned  to  each  node  j  G  AT.  In  this  illustration,  we  use  the  continuous  ex¬ 
ponential  distribution  with  rate  parameter  p  and  cumulative  distribution  function 
(c.d.f) 

{1  —  exp(— px)  x  >  0, 

0  x  <  0 

The  rate  parameter  for  node  j,  pj ,  is  recorded  in  Table  4.5. 

Table  4.5  Rate  parameters  chosen  for  the  exponential  distributions  used  for  example  2. 


Node 

Pj 

1 

6666.67 

2 

1538.46 

3 

1322.75 

4 

1574.80 

5 

1754.39 

6 

10000.00 

7 

4347.83 

8 

2222.22 

9 

909.09 

10 

7692.31 

The  expected  value  of  the  random  vapor  concentration,  E[Vj]  for  each  node 
j  G  AT  is  computed  by 

m  1  =  1/n-  (4.5) 

These  expected  values  are  used  to  determine  the  reward  value  r3  for  each  node.  Both 
of  these  values  are  recorded  in  Table  4.6. 
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Table  4.6  Expected  vapor  concentrations  and  rewards  for  example  2  (a*  =  6.0  x  10  4). 


Node 

£K](*io-4) 

ri 

1 

1.50 

1 

2 

6.50 

0 

3 

7.56 

0 

4 

6.35 

0 

5 

5.70 

1 

6 

1.00 

1 

7 

2.30 

1 

8 

4.50 

1 

9 

11.00 

0 

10 

1.30 

1 

Now  that  the  expected  value  and  expected  reward  for  searching  each  node 
has  been  calculated,  the  same  MATLAB®  code  used  in  section  4.2  is  implemented 
to  find  the  optimal  path  and  total  reward  for  sampling  the  areas.  Again,  since  the 
algorithm  seeks  to  find  the  path  of  shortest  time  throughout  the  site,  the  travel  times 
have  not  changed,  and  the  same  starting  node  s  =  5  is  chosen.  The  optimal  path  is 
identical  to  the  optimal  path  identified  in  the  case  of  static  and  deterministic  vapor 
concentrations.  Thus,  the  solution  is  again 


^  =  {5,8,10,7,4,1,2,3,6,9} 


with  a  total  search  time  of  62.85  minutes.  The  solution  is  illustrated  in  Figure  4.4 
where  the  bold  nodes  indicate  locations  containing  vapor  concentration  levels  below 
the  threshold. 

It  is  important  to  note  the  total  reward  has  not  changed  from  6  but  the  six 
areas  are  not  identical.  This  is  due  to  the  stochastic  behavior  of  the  secondary  vapor 
concentrations  at  each  node.  Since  the  total  reward  in  this  case  is  based  on  expected 
values  of  the  vapor  concentrations,  searching  the  site  may,  in  fact,  yield  an  actual 
reward  value  that  differs  from  this  reward.  In  this  example,  the  areas  in  which  a 
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6 


9 


The  next  section  contains  an  illustration  of  the  case  with  dynamic  and  deter¬ 
ministic  vapor  concentrations.  Unlike  the  two  previous  examples,  the  vapor  con¬ 
centrations  evolve  over  time,  forcing  the  search  crew  to  re-assess  their  route  in  a 
time-adaptive  search. 


4-4  Example  3:  Dynamic  and  Deterministic  Vapor  Concentrations 

Just  as  in  the  case  of  section  4.2,  we  assume  vapor  concentrations  can  be 
calculated  at  initial  time  t0  using  a  deterministic  formula.  Using  the  same  site  and 
travel  times  from  the  previous  two  examples,  we  again  use  the  advection-diffusion 
Equation  (3.1)  discussed  in  Chapter  3  to  determine  vapor  concentrations  for  each 
area  of  the  site. 
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Consider  again  the  ten-node  site  depicted  in  Figure  4.2.  Let  the  search  crew 
begin  at  a  node  chosen  arbitrarily,  say  node  5  (i.e.,  the  starting  node  s  —  5).  Using 
Equation  (3.1),  the  vapor  concentration  for  the  starting  node  v5(t0)  is  computed 
using  the  initial  time  to  =  7  min  (i.e.,  time  since  the  chemical  agent  release). 

Using  the  algorithm  outlined  in  Chapter  3,  subsection  3.2.3,  and  the  time 
matrix  T.  the  initialization  step  yields  the  results  in  Table  4.7 


Table  4.7  Results  of  initialization  step  for  dynamic/deterministic  example. 


Node 

t  (sec) 

Vj{t) 

rj(i) 

5 

420 

0.00065 

0 

The  set  77  remains  empty  and  the  set  0  =  {5}.  In  step  2,  the  vapor  concen¬ 
tration  at  nodes  adjacent  to  node  5  are  examined  using  the  current  time  t0  plus  the 
travel  time  t,h3  from  node  5  to  those  adjacent  nodes  j  G  A f%  (t  —  to  +  Uj).  Recall,  the 
time  tij  includes  5  minutes  required  for  obtaining  an  instrument  reading  at  node  j. 
An  indicator  value  of  Tj(t0)  =  1  is  assigned  to  node  j  if  v3(t)  <  v*.  The  results  of 
iteration  1  are  displayed  in  Table  4.8.  The  anticipated  reward  of  proceeding  to  each 
node  adjacent  to  node  5  is  0  so  the  arc  of  shortest  time  is  chosen  (i.e.,  the  arc  from 
5  to  8).  The  following  updates  are  recorded: 

t  =  420  +  369.6  =  789.6; 

77  =  0; 

0  =  0U{8}  =  {5,8}. 

Now  the  vapor  concentration  for  each  node  adjacent  to  node  8  is  calculated 
for  t  =  789.6  +  t8j  for  all  j  G  Af8.  Though  the  adjacent  set  Af8  consists  of  the 
nodes  {5,  6,  9, 10},  only  nodes  6,  9,  and  10  are  considered  since  node  5  was  previously 
sampled.  That  is,  the  crew  cannot  return  to  the  node  from  which  it  has  just  departed. 
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Table  4.8  Results  of  iteration  1  for  dynamic/deterministic  example  (v*  =  6.0  x  10  4). 


Table  4.9 


Node 

t  (sec) 

Vj{t) 

ri 

1 

858.6 

0.0036 

0 

2 

822.6 

0.0035 

0 

3 

832.8 

0.0036 

0 

7 

825.6 

0.0064 

0 

8 

789.6 

0.0054 

0 

•ation  2  for  dynamic/deterministic 

Node 

t  (sec) 

Vj(t) 

rj(t) 

6 

1220.4 

0.0048 

0 

9 

1266.0 

0.0080 

0 

10 

1123.8 

0.0070 

0 

-1-4  \ 


After  iteration  2,  node  10  is  chosen  as  the  next  node  for  the  crew  to  search 
because  the  reward  for  all  nodes  adjacent  to  node  8  is  0  and  the  least  time  path 
is  from  8  to  10.  The  nodes  A/j0  =  {4,  7,  8}  are  examined  next.  However,  node  8 
is  momentarily  excluded  from  consideration  since  it  has  just  been  searched.  The 
algorithm  continues  for  12  iterations  and  terminates  when  the  allotted  time  of  90 
minutes  (5400  seconds)  is  reached.  The  results  of  each  iteration  are  recorded  in 
Table  4.10  where  the  node  in  bold  is  the  next  node  to  be  searched. 
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Table  4.10 


Iterations  3-12  for  dynamic/deterministic  example  (a*  =  6.0  x  10  4). 


Node 

t  (sec) 

Vj(t) 

rj(t) 

Iteration  3 

4 

1671.6 

0.0035 

0 

7 

1569.6 

0.0056 

0 

Iteration  4 

4 

1989.6 

0.0024 

0 

5 

1975.2 

0.0028 

0 

Iteration  5 

1 

2413.8 

0.00092 

0 

2 

2377.8 

0.0010 

0 

3 

2388.0 

0.0012 

0 

8 

2344.8 

0.0022 

0 

Iteration  6 

6 

2775.6 

0.00091 

0 

9 

2821.2 

0.0013 

0 

10 

2679.0 

0.0014 

0 

Iteration  7 

4 

3226.8 

0.00041 

1 

7 

3124.8 

0.00072 

0 

Iteration  8 

1 

3618 

0.00016 

1 

3 

3720.8 

0.00017 

1 

7 

3646.8 

0.00033 

1 

Iteration  9 

5 

4052.4 

0.00014 

1 

10 

4092.6 

0.00019 

1 

Iteration  10 

4 

4640.4 

0.000050 

1 

8 

4426.8 

0.00011 

1 

Iteration  11 

5 

4796.4 

0.000047 

1 

6 

4857.6 

0.000044 

1 

9 

4903.2 

0.000059 

1 

Iteration  12 

6 

5506.2 

0.000017 

1 
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The  algorithm  terminates  after  12  iterations,  yielding  the  path 


^  =  {5,  8, 10,  7, 5,  8, 10, 4,  7, 10, 8,  9} 

and  reduced  set 

77  =  {4,  7,8,9,10}. 

Area  6  would  be  the  next  area  to  search;  however,  the  time  required  to  travel  to  6  and 
take  the  measurement  exceeds  the  allotted  time  by  1.77  minutes  so  the  algorithm 
terminates  at  node  9.  The  total  reward  for  this  example  is  5,  where  areas  not 
requiring  protective  gear  are  the  nodes  j  E  1Z.  The  total  time  required  to  complete 
the  search  (including  obtaining  instrument  readings)  is  81.72  minutes.  Figure  4.5 
graphically  displays  the  search  policy  with  bold  nodes  denoting  those  areas  not 
requiring  protective  gear.  The  numbers  along  each  arc  indicate  the  order  in  which 
the  arcs  are  added  to  the  solution  (i.e.,  the  order  of  the  sampling). 


Figure  4.5  Optimal  sampling  path  when  vapor  concentrations  are  dy¬ 
namic  and  deterministic. 
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It  is  important  to  note  the  state  of  the  system  at  the  terminating  time.  Using 
Equation  (3.1),  the  vapor  concentration  and  reward  value  at  each  node  is  calculated 
for  t  =  81.72  min  and  these  values  are  tabulated  in  Table  4.11.  According  to  this 
information,  all  10  areas  of  the  site  can  safely  operate  without  the  protective  gear. 
However,  the  search  crew  cannot  reach  all  of  those  areas  after  the  time  when  the 
vapor  concentration  decreases  below  the  threshold  and  terminate  the  search  within 
the  allotted  time  of  90  minutes. 


Table  4.11  Vapor  concentration  and  reward  for  each  node  at  algorithm’s  termination 
(v*  =  6.0  x  10-4). 


Node 

Uj(81.72) 

rj(t) 

1 

0.000023 

1 

2 

0.000025 

1 

3 

0.000030 

1 

4 

0.000034 

1 

5 

0.000040 

1 

6 

0.000042 

1 

7 

0.000052 

1 

8 

0.000054 

1 

9 

0.000059 

1 

10 

0.000057 

1 

In  the  following  section,  we  provide  an  example  for  the  case  of  dynamic  and 
stochastic  vapor  concentration  levels.  We  employ  the  algorithm  of  Chapter  3,  sub¬ 
section  3.2.4  to  identify  a  policy  for  directing  a  search  crew  through  the  site,  by 
searching  those  areas  most  likely  to  have  a  vapor  concentration  below  the  threshold 
level. 

4-5  Example  4:  Dynamic  and  Stochastic  Vapor  Concentrations 

Finally,  we  provide  an  illustration  of  the  case  of  dynamic  and  stochastic  va¬ 
por  concentrations.  More  representative  of  the  real-world  problem,  we  allow  the 
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secondary  vapor  concentrations  to  evolve  over  time  according  to  a  time-dependent 
probability  distribution. 

Consider  the  notional  ten-node  site  depicted  in  Figure  4.2.  Again,  the  time  to 
travel  from  node  i  to  node  j  remains  hxed  and  the  matrix  of  times  T  used  for  the 
previous  examples  is  used  here.  Continuous  probability  distributions  and  parameters 
are  assigned  to  each  node  j  G  A/".  For  this  example,  we  use  the  same  distributions 
and  initial  parameters  used  in  example  2  of  section  4.3.  Recall,  we  assume  each  node 
is  governed  by  the  exponential  distribution  with  their  initial  rate  parameters  listed 
in  Table  4.5. 

In  the  case  of  dynamic  and  stochastic  secondary  vapor  concentrations,  the 
probability  distributions  change  over  time.  Thus,  for  illustrative  purposes,  we  choose 
the  relationship 

Hj{t  +  tij)  —  (t  +  tij)  x  x  0.01  (4.6) 

to  model  the  evolution  of  the  distribution  over  time.  Thus,  at  each  iteration  of  the 
algorithm,  the  parameters  are  updated  according  to  the  time  at  which  the  search 
crew  arrives  at  node  j  G  A/". 

The  algorithm  presented  in  subsection  3.2.4  of  Chapter  3  is  implemented  to 
find  the  path  for  this  case.  As  opposed  to  the  algorithm  employed  in  the  previous 
section,  it  is  the  probability  that  a  future  node  has  a  vapor  concentration  below  the 
threshold  that  drives  the  solution  identified  at  the  termination  of  this  algorithm. 
We  again  choose  s  =  5  as  the  starting  location  at  time  t0  =  420  seconds.  Using 
the  probability  distribution  and  parameter  for  node  5,  a  realization  1^(420)  of  the 
random  variable  Vs(420)  is  obtained  by  drawing  a  number  from  the  appropriate 
exponential  population  (i.e.,  14(420)  ~  Exp(p5( 420))). 

Table  4.12  displays  the  results  of  the  initialization  step.  At  time  t0  =  420  sec, 
there  is  a  0.65  probability  that  the  vapor  concentration  at  node  5  is  less  than  the 
threshold  v*  =  0.0006  mg/m3.  The  realization  of  the  vapor  concentration  at  node 
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5  at  time  420  sec  is  0.0011  mg/m3.  The  realization  value  represents  the  instrument 
reading  that  will  be  taken  at  that  time  t  in  area  5  during  the  implementation  of  the 
path. 


Table  4.12  Initialization  step  for  dynamic/stochastic  example. 

Node  t  (sec)  P(Vj(t)  <  v*) 

5  420.0  1754.39  (765 


A  reward  value  of  0  is  assigned  to  node  5  for  the  current  time  and  all  nodes 
adjacent  to  node  5  are  considered  in  the  next  step.  The  probability  that  the  va¬ 
por  concentration  is  below  the  fixed  threshold  for  all  nodes  adjacent  to  node  5  is 
calculated.  Based  on  these  probabilities,  recorded  in  Table  4.13,  the  next  node  is 
determined. 


Table  4.13  Results  of  step  1  for  dynamic/stochastic  example. 


Node 

t  (sec) 

wit)  pm)  < v*) 

1 

858.6 

775.19 

0.372 

2 

822.6 

186.92 

0.106 

3 

832.8 

158.73 

0.091 

7 

825.6 

526.32 

0.271 

8 

789.6 

281.69 

0.155 

The  node  chosen  to  be  searched  from  node  5  is  node  1  since  it  has  the  highest 
probability  of  having  a  vapor  concentration  below  the  threshold.  Continuing  in  this 
fashion  the  algorithm  terminates  after  twelve  iterations.  The  resulting  path  is 

^  =  {5,1,5,7,10,8,6,3,5,1,4, 10} 

and  this  search  is  completed  in  88.61  minutes. 

Drawing  random  numbers  from  the  family  of  exponential  distributions  for  each 
node  at  each  time,  we  determine  the  vapor  concentration  Vj(t)  for  all  j  G  at  the 
time  they  are  searched.  The  values  are  recorded  in  Table  4.14. 
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Table  4.14 


Vapor  concentrations  for  nodes  in  ip  for  dynamic/stochastic  example  (u*  = 
6.0  x  10"4). 


Node 

t  (sec) 

Vj(t) 

rj(t) 

5 

420 

0.0011 

0 

1 

858.6 

0.000012 

1 

5 

1297.2 

0.000106 

1 

7 

1702.8 

0.000172 

1 

10 

2148.6 

0.000103 

1 

8 

2482.8 

0.000113 

1 

6 

2913.6 

0.000172 

1 

3 

3345.6 

0.000179 

1 

5 

3758.4 

0.000609 

0 

1 

4197.0 

0.000693 

0 

4 

4588.2 

0.000803 

0 

10 

5136.0 

0.000001 

1 

Note,  due  to  the  dynamic  and  stochastic  behavior  of  the  secondary  vapor  con¬ 
centrations,  nodes  that  initially  have  concentrations  below  the  threshold  may  become 
contaminated  at  later  times.  The  total  reward  of  this  policy  at  the  termination  of 
the  allotted  search  time  is  5.  Though  the  sum  of  the  reward  values  in  Table  4.14 
is  8,  one  of  those  nodes  (node  10)  is  repeated  and  two  (nodes  1  and  5)  have  vapor 
concentrations  that  increase  above  the  threshold  v*  when  searched  again.  Hence, 
the  reduced  set  at  the  termination  is 


77  =  {3,6,7,8,10}. 

Thus,  these  5  areas  are  the  areas  most  likely  to  operate  safely  without  protective 
gear. 

Figure  4.6  is  the  graphical  representation  of  the  path  tp,  where  the  bold  nodes 
represent  areas  where  protective  gear  is  not  required.  The  numbers  along  each  arc 
indicate  the  order  in  which  the  arc  is  added  to  the  solution. 
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Figure  4.6  Optimal  sampling  path  when  vapor  concentrations  are  dy¬ 
namic  and  stochastic. 

4-6  Comparison  of  Static  and  Dynamic  Deterministic  Solutions 

The  examples  of  sections  4.2  and  4.4  provide  the  basis  for  analyzing  the  im¬ 
pact  of  ignoring  (or  incorporating)  time  into  the  mathematical  models  developed  in 
Chapter  3.  We  show  there  is  a  difference  between  the  case  where  vapor  concentration 
levels  remain  static  and  when  they  evolve  over  time. 

Table  4.15  displays  the  solution  to  example  1  (static  and  deterministic)  and 
example  3  (dynamic  and  deterministic).  The  optimal  path  for  example  1 

^  =  {5,8,10,7,4,1,2,3,6,9} 

directs  the  search  crew  through  all  ten  nodes,  reaching  each  one  exactly  once,  in 
t*  =  62.85  minutes.  The  total  reward,  found  from  deterministically  calculating  the 
vapor  concentration  for  each  node  j  E  AT,  is  6.  They  should  find  that  nodes  in  the 
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reduced  set 


77  =  {1,2,3,6,8,10} 

have  vapor  concentrations  below  the  fixed  threshold. 

The  optimal  path  of  example  3, 

i\)  =  {5,  8, 10,  7, 5,  8, 10, 4,  7, 10, 8,  9} 

requires  r*  =  81.72  minutes  to  complete,  and  never  reaches  nodes  1,  2,  3,  or  6. 
The  neglect  of  these  nodes  is  due  to  the  dynamic  nature  of  the  secondary  vapor 
concentrations.  In  this  case,  the  concentration  of  the  node  at  the  crew’s  arrival  time 
drives  the  algorithm.  For  all  nodes  with  a  concentration  below  the  threshold,  the 
algorithm  always  chooses  the  vapor  concentration  closest  to  the  threshold.  If  none 
of  the  nodes  will  have  a  concentration  below  the  threshold,  the  algorithm  chooses 
the  closest  (in  distance  or  time)  node.  This  also  allows  nodes  to  be  searched  if  the 
vapor  concentration  decreases,  even  though  it  has  been  previously  searched. 

Note,  in  example  1  (static/deterministic  vapor  concentrations),  27.15  minutes 
remain  in  the  time  allotted  for  the  search.  We  know  vapor  concentrations  evolve 
over  time,  so  those  extra  minutes  are  valuable  in  allowing  the  search  crew  to  reach 
more  nodes  at  later  times  to  record  any  change  in  concentration,  either  positive  or 
negative. 

Figures  4.7  and  4.8  show  the  temporal  evolution  of  the  vapor  concentrations 
for  each  node  throughout  the  duration  of  the  search  for  example  3.  The  vapor 
concentration  from  the  static  case  in  example  1  can  also  be  seen  on  the  plot.  They 
are  the  vapor  concentrations  at  time  420  seconds.  It  is  clear  that  time  directly 
impacts  the  vapor  concentrations  and  thus,  must  be  considered  when  developing  a 
search  policy  that  will  optimally  route  a  crew. 

In  example  3,  only  fifty  percent  of  the  nodes  were  found  to  have  a  vapor 
concentration  less  than  the  threshold,  whereas  in  example  1,  sixty  percent  of  the 
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Table  4.15  Comparison  of  solutions  to  the  static/deterministic  and  dy¬ 
namic  /deterministic  examples. 


Static 

Dynamic 

Node 

ri 

Node 

rAT) 

5 

0 

5 

0 

8 

1 

8 

0 

10 

1 

10 

0 

7 

0 

7 

0 

4 

0 

5 

0 

1 

1 

8 

0 

2 

1 

10 

0 

3 

1 

4 

1 

6 

1 

7 

1 

9 

0 

10 

1 

8 

1 

9 

1 

Total  Tune  (min) 

r* 

Total  Time  (min) 

r*(r*) 

62.85 

6 

81.72 

5 

nodes  no  longer  required  protective  gear.  However,  only  two  of  those  nodes  have 
vapor  concentrations  below  the  threshold  in  the  example  solutions.  Thus,  had  the 
effect  of  time  been  ignored,  personnel  in  areas  1,  2,  3,  and  6  may  have  been  harmed 
by  not  donning  their  protective  equipment.  The  safe  nodes  of  example  3  (nodes  4, 
7,  8,  9,  and  10)  may  be  more  accurate  since  they  account  for  the  evolution  of  the 
vapor  concentrations  over  time. 

In  the  next  section,  we  compare  the  results  of  examples  2  and  4.  This  allows  us 
to  examine  the  impact  of  temporal  evolution  when  we  assume  vapor  concentrations 
are  stochastic  rather  than  deterministic. 


f.,7  Comparison  of  Static  and  Dynamic  Stochastic  Solutions 

The  comparison  of  examples  2  and  4  provide  an  interesting  analysis  of  the  effect 
of  time  of  the  secondary  vapor  concentrations  when  they  are  assumed  to  behave 
according  to  continuous  probability  distributions.  In  example  2,  the  probability 
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Figure  4.7  Evolution  of  vapor  concentrations  over  time  from  example 
3,  nodes  1  through  5. 


distributions  are  static  and  in  example  4,  they  are  time-variant.  Recall,  the  result 
of  example  2  yields  a  path 


^  =  {5,8,10,7,4,1,2,3,6,9} 


directing  the  search  crew  to  all  ten  areas,  reaching  each  one  exactly  once.  In  example 
4,  when  time  is  incorporated  into  the  methodology,  the  resulting  path 


^  =  {5,1,5,7,10,8,6,3,5,1,4, 10} 


directs  the  search  crew  to  just  8  of  the  10  nodes.  However,  in  example  4,  the  search 


requires  88.61  minutes  as  opposed  to  the  62.85  required  for  example  2.  The  longer 
search  time  in  the  dynamic  case  is  a  result  of  allowing  previously  searched  nodes  to 
be  searched  again.  Since  the  expected  vapor  concentrations  are  time-dependent  there 
exists  a  nonzero  probability  that  the  vapor  concentration  at  a  node  may  increase  (or 
decrease)  over  time.  Thus,  several  of  the  nodes  are  searched  more  than  once  while 
others  are  never  searched  at  all. 

In  both  examples,  the  vapor  concentrations  behave  according  to  exponential 
probability  distributions.  It  is  the  probability  that  a  vapor  concentration  is  below 
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Figure  4.8  Evolution  of  vapor  concentrations  over  time  from  example 
3,  nodes  6  through  10. 

the  threshold  that  drives  the  node  selection  in  the  algorithm  in  example  4.  Thus,  it 
directs  the  crew  as  follows: 

i\)  =  {5,1,5,7,10,8,6,3,5,1,4,10}. 

The  reward  for  case  4  stems  from  searching  nodes  in  the  reduced  set 

77  =  {3,6,7,8,10}. 

Whereas  the  solution  to  example  2  is  to  search  the  nodes  in  the  following  order: 

^  =  {5,8,10,7,4,1,2,3,6,9} 
with  a  reward  for  searching  nodes 

77  =  {1,5,6,7,10}. 
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Table  4.16  Comparison  of  solutions  to  the  static/stochastic  and  dynamic/stochastic 
examples. 


Static 

Dynamic 

Node 

r3 

Node 

rj(T) 

5 

1 

5 

0 

8 

1 

1 

0 

10 

1 

5 

0 

7 

1 

7 

1 

4 

0 

10 

1 

1 

1 

8 

1 

2 

0 

6 

1 

3 

0 

3 

1 

6 

1 

5 

0 

9 

0 

1 

0 

4 

0 

10 

1 

Total  Time  (min)  r* 

Total  Time  (min) 

r*(r*) 

62.85 

6 

88.61 

5 

The  case  assuming  stochastic  and  dynamic  vapor  concentrations  suggests  that 
50  percent  of  the  nodes  no  longer  require  protective  gear.  If  the  policy  for  example 
2  had  been  followed,  they  would  believe  that  60  percent  of  the  areas  can  operate 
safely  without  protective  gear;  thus  jeopardizing  the  personnel  in  those  areas  where 
the  progression  of  time  may  cause  the  vapor  concentrations  to  contaminate  those 
areas  that  may  have  previously  been  deemed  “safe.” 

4-8  Summary  of  Results 

The  four  numerical  examples  provided  in  this  chapter  demonstrate  the  impact 
of  increasing  the  complexity  of  the  problem  on  the  solution.  In  the  first  and  most 
simplistic  case,  we  determine  that  the  optimal  path  is  one  that  samples  every  node 
exactly  once  in  the  minimum  amount  of  time.  However,  this  case  does  not  capture 
the  true  behavior  of  secondary  vapor  concentrations  that  are  known  to  evolve  over 
time  and  are  subject  to  stochastic  elements  (e.g.,  wind,  temperature,  etc.).  In  the 
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real  problem,  it  may  be  beneficial  to  re-examine  an  area  if  it  is  likely  the  vapor 
concentration  has  changed.  Thus,  the  state  of  the  system  in  the  first  example  may 
not  accurately  represent  the  true  state,  resulting  in  a  suboptimal  solution  of  the 
real  problem.  The  second  case  assumes  vapor  concentrations  are  governed  by  time- 
invariant  probability  distributions.  The  optimal  path  is  identical  to  the  path  of  the 
previous  example,  but  the  total  reward  is  different  due  to  the  stochastic  nature. 

The  last  two  numerical  illustrations  demonstrate  the  impact  time  has  on  the 
solution.  In  example  3,  we  assume  the  vapor  concentrations  can  be  computed  by  an 
equation  that  depends  on  time.  The  crew  does  not  search  every  node  in  the  optimal 
solution  and  searches  some  more  than  once.  The  solution  procedure  accounts  for 
the  evolution  of  the  vapor  concentrations  over  time  and  chooses  where  to  direct  the 
crew  based  on  the  calculated  concentration  at  their  time  of  arrival.  Thus,  the  state 
of  the  system  in  example  3  is  more  representative  of  reality  than  either  example  1 
or  2.  Finally,  example  4  attempts  to  capture  the  true  nature  of  the  problem.  The 
vapor  concentration  at  each  node  evolves  according  to  a  time-dependent  probability 
distribution.  As  time  progresses,  the  probability  a  node  has  a  vapor  concentration 
below  (or  above)  the  threshold  also  evolves.  Thus,  in  these  illustrations,  the  resulting 
path  accounts  for  the  effect  time  has  on  the  chemical  agent  as  well  as  the  effect  time 
has  on  the  probability  distribution. 

Consider  the  result  of  example  1,  the  case  of  static  and  deterministic  vapor 
concentrations.  The  optimal  path  is  one  that  reaches  each  node  exactly  once.  Com¬ 
paring  this  result  with  that  of  the  dynamic  and  deterministic  example  emphasizes 
the  impact  of  incorporating  temporal  evolution  into  the  model.  In  example  3,  the 
optimal  path  directs  the  search  crew  to  some  nodes  more  than  once  and  others  not 
at  all.  Due  to  the  temporal  evolution  of  vapor  concentrations,  some  areas  are  more 
critical  than  others  and  require  being  searched  several  times,  and  those  that  will  not 
change  enough  to  warrant  a  search  are  excluded  from  the  solution.  Similar  compar¬ 
isons  can  be  made  with  examples  2  and  4,  the  cases  considering  stochastic  vapor 
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concentrations.  Observing  the  results  of  both  dynamic  cases  and  comparing  them 
with  their  static  counterparts  demonstrates  that  ignoring  the  temporal  evolution 
of  the  vapor  concentrations  could  lead  to  suboptimal  results  that  overestimate  (or 
underestimate)  the  number  of  areas  operating  safely  without  protective  gear.  For 
example,  if  the  optimal  path  of  example  1  is  followed,  four  areas  would  be  operating 
dangerously  without  protective  gear  because  in  example  3,  those  four  areas  are  found 
to  be  hazardous  at  later  times.  Also,  three  of  the  areas  would  be  inhibited  by  the 
cumbersome  gear,  when  they  could  be  operating  safely  without  it.  Likewise,  the 
stochastic  cases  can  be  compared  with  a  similar  outcome. 

It  is  important  to  note  the  examples  in  this  chapter  provide  numerical  illustra¬ 
tions  based  on  several  restrictive  assumptions.  A  possible  explanation  for  the  results 
observed  in  these  illustrations  is  that  we  assumed  independent  vapor  concentration 
samples  at  each  area  over  time.  This  may  explain  why  we  observed  higher  vapor 
concentrations  in  some  areas  at  a  point  in  time  when  it  was  previously  below  the 
threshold.  In  order  to  relax  this  assumption,  we  need  a  stochastic  process  in  time  and 
space  describing  the  vapor  concentration  evolution.  In  the  next  chapter  we  provide 
a  summary  of  this  thesis,  concluding  remarks,  and  propose  areas  of  future  work. 
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5.  Conclusions  and  Future  Research 


The  threat  of  terrorist  attacks  with  biological  or  chemical  weapons  has  become 
a  serious  consideration  for  nations  and  their  military  forces.  Modern  research  is  un¬ 
derway  for  identifying  emergency  response  procedures  for  specific  biological  weapons 
and  for  the  incorporation  of  new  sensor  technologies  and  instruments  for  identifying 
and  quantifying  chemical  contamination.  However,  methods  for  effectively  sampling 
(in  an  optimal  manner)  chemical  hazard  areas  subsequent  to  a  chemical  weapons 
attack  lag  far  behind.  This  thesis  has  provided  a  methodological  framework  for 
incorporating  real-time  information  for  the  purpose  of  optimally  sampling  a  site 
contaminated  by  an  initial  chemical  attack  and  hazardous  secondary  chemical  vapor 
concentrations. 

Three  distinct  cases  of  the  problem  were  considered  before  the  fourth  (and 
most  realistic  case)  was  examined.  These  cases  were:  static  and  deterministic  vapor 
concentration  levels,  static  and  stochastic  vapor  concentration  levels,  dynamic  and 
deterministic  vapor  concentration  levels,  and  dynamic  and  stochastic  vapor  con¬ 
centration  levels.  For  each  of  the  four  cases,  the  objective  function  and  relevant 
constraints  were  defined  and  a  solution  methodology  presented.  The  first  two  cases 
ignored  the  temporal  evolution  of  vapor  concentrations  and  focused  on  identifying 
the  optimal  path  through  a  static  system  to  maximize  the  reward  (i.e.,  the  number 
of  areas  no  longer  requiring  protective  gear).  The  final  two  cases  incorporated  tem¬ 
poral  evolution  of  the  chemical  vapor  concentrations  to  more  accurately  model  the 
real  world  problem.  The  results  of  this  thesis  represent  a  significant  contribution  to 
the  military  operations  research  and  bio-chemical  defense  communities  by  providing 
a  method  that  can  be  implemented  quickly  and  efficiently  in  response  to  a  chemical 
attack  on  a  military  installation. 

The  methodology  for  each  case  was  illustrated  on  four  numerical  example  prob¬ 
lems  to  identify  the  optimal  path  that  results  in  the  maximum  reward  within  the 
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allotted  time  for  the  search.  The  parameters  were  defined,  and  the  algorithm  de¬ 
veloped  for  each  specific  case  was  implemented  to  obtain  those  results.  Finally,  the 
result  of  the  two  static  examples  were  compared  with  their  respective  dynamic  coun¬ 
terparts  to  highlight  the  impact  of  the  temporal  evolution  of  vapor  concentrations 
on  the  solution.  It  is  clear  that  incorporating  this  evolution  into  the  model  resulted 
in  solutions  that  captured  more  of  the  real  problem’s  dynamics,  thus  maximizing  the 
number  of  areas  in  an  attacked  site  that  can  operate  safely  without  the  cumbersome 
protective  gear.  Ignoring  the  temporal  evolution  of  the  vapor  concentrations  may 
result  in  removing  the  protective  gear  requirement  in  areas  that  may  subsequently 
become  contaminated,  thus  exposing  humans  to  an  unnecessary  hazard. 

To  realistically  implement  the  techniques  presented  in  this  thesis,  all  of  the 
parameters  relevant  to  the  methodology  must  be  known.  This  thesis  covers  a  small 
part  of  a  much  larger  problem:  estimating  future  downwind  hazard  areas  when  only 
detector /sensor  data  (e.g.,  human  input,  M-8/M-9  paper,  etc.)  is  known  and  using 
those  inputs  to  also  estimate  the  source  location  and  strength  of  the  chemical  agent 
release.  The  last  part  of  the  problem  is  to  estimate  the  area  of  contamination.  We 
assume  this  information  is  known  a  priori ,  allowing  us  to  optimally  direct  a  search 
crew  through  the  site  to  sample  areas  likely  to  have  vapor  concentrations  below  a 
fixed  threshold,  indicating  the  area  is  safe. 

Despite  the  fact  that  we  assume  the  dynamics  of  the  problem,  the  results  of 
this  thesis  provide  insight  into  data  requirements  for  optimally  sampling  a  chemical 
hazard  area.  For  example,  a  model  for  estimating  the  spatiotemporal  evolution  of 
secondary  vapor  concentrations  is  critical  to  providing  information  regarding  the 
concentration’s  probability  distribution  for  each  area.  Employing  a  realistic  model 
provides  a  means  to  identify  areas  likely  (or  unlikely)  to  be  contaminated.  These 
distributions  may  subsequently  be  applied  to  the  stochastic  cases  of  sections  3.2.2 
and  3.2.4.  Once  the  information  is  obtained,  one  can  identify  the  areas  that  can 
safely  operate  without  protective  gear. 
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Another  restrictive  assumption  of  this  thesis  is  the  use  of  static  and  determinis¬ 
tic  travel  times  between  adjacent  areas  of  the  site.  Traffic,  obstacles,  and  other  such 
impediments  are  likely  to  delay  the  search  crew  or  prohibit  it  from  taking  a  path 
altogether.  Incorporating  dynamic  and  stochastic  travel  times  will  result  in  solutions 
with  a  least  expected  time  path  with  a  maximum  reward.  Relaxing  this  assump¬ 
tion  also  leads  to  possible  research  in  incorporating  real-time  information  into  the 
methodology.  For  example,  as  routes  open  or  close  clue  to  unforseen  circumstances, 
the  solution  path  can  be  updated  to  accommodate  these  changes. 

Relaxing  the  assumption  of  a  single  search  crew  would  also  result  in  different 
methodology.  If  more  than  one  crew  were  available,  a  higher  number  of  areas  may 
be  searched  within  the  allotted  time,  resulting  in  more  areas  operating  in  the  correct 
level  of  protective  gear.  Having  a  number  of  crews  equal  to  the  number  of  areas 
would  be  ideal  because  the  problem  would  reduce  to  finding  the  shortest  path  from 
each  crew’s  location  to  a  specific  area. 

Finally,  another  area  for  future  research  involves  incorporating  real-time  data 
into  the  methodology  to  better  model  the  behavior  of  the  secondary  vapor  concen¬ 
trations.  This  is  due  to  the  temporal  and  spatial  evolution  as  well  as  stochastic 
elements  (e.g.,  wind  direction  and  velocity,  temperature,  humidity,  etc.)  affecting 
the  behavior.  Once  future  downwind  hazard  areas  can  be  identified  and  the  inputs 
required  for  estimating  parameters  for  those  areas  can  be  determined,  they  can  be 
incorporated  into  the  techniques  developed  in  this  thesis  providing  an  optimal  path. 

The  main  results  of  this  thesis  reach  far  beyond  military  applications  only. 
Hazardous  vapor  concentrations  resulting  from  industrial  accidents  may  contaminate 
surrounding  communities,  threatening  the  lives  of  civilians  within  some  radius  of  the 
contamination.  Having  a  procedure  to  optimally  sample  the  affected  area  as  quickly 
as  possible  (to  identify  which  residents  should  be  evacuated)  is  invaluable  to  saving 
lives.  Also,  should  a  chemical  agent  be  released  over  a  city  or  region,  such  a  procedure 
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may  directly  help  emergency  crews  identify  contaminated  areas  and  evacuate  them 
in  an  efficient  and  timely  manner. 
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Appendix  A.  Static/Deterministic  Code 


7:/.y.y//:/.7.y:/:/.y.7.y:/:/.y.7:/:/.7.y:/:/.y.y:/:/:/.y.7:/:/.7.y:/:/.y.y:/;/.7.y.y:/;/.7.y:/.7.y.y:/:/.y.7.,/:/.7.y.y:/:/.7.,/:/.7.y.y:/:/. 


7.  Computing  STATIC  vapor  concentrations  for  (x,y)  coordinates  7, 
7.  7. 
7,  The  purpose  of  this  MATLAB  code  is  to  compute  the  static  and  7« 
7»  deterministic  vapor  concentrations  for  the  network  used  in  7« 
7.  this  thesis.  7« 


7. 


7. 


7.  Author:  Jennifer  R.  Plourde,  AFIT/ENS  G0R-05M  7« 
7.  Date:  1  Feb  05  7« 
7,  Last  revised:  1  Feb  05  7. 
7.  References:  Kathirgamanathan,  P.,  McKibbin,  R.  and  7« 
7,  R.I.  McLachlan  (2003).  Source  release-rate  estimation  7« 
7.  of  atmospheric  pollution  from  a  non-steady  point  source  7« 
7,  -  Part  1:  Source  ate  a  known  location.  Res.  Lett.  Inf.  % 
7,  Math.  Sci.  vol  5,  71-84.  % 
7.  7. 


7o7.7o7o7o7o7.7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7.7o7o7o7o7o7o7o707o7.7o7o7o7o7o7o7o7«7o7o7o7o7.7o7o7o7o7o7o7o7o7.707o7o7o7o7o7o7.7o7o7.7o7o707o 


7o7.7o7o7o7o7o7o7o7o7o7o7o7o7o7o7.7o7o7.7o7.7o  variable  definitions  7o7o7o7o7o7o7o7o7o7o7o7o7.7o7o7o7o7.7o7o7o7o7o7o7. 


7,  Inputs  required:  7. 
7.  nel:  number  of  nodes  7. 
7.  tm:  (nel  x  nel)matrix  of  feasible  travel  times  between  nodes  7. 
7.  X:  (nel  x  1)  vector  of  x-coordinates  7. 
7,  Y:  (nel  x  1)  vector  of  y-coordinates  7. 
7.  Z:  (nel  x  1)  vector  of  z-coordinates  (z  =  0,  for  this  thesis)  7. 
7,  V :  vapor  threshold  value  7. 
7.  Q :  total  mass  release  (kg)  7. 
7.  K_x,  K_y,  K_z:  eddy  dif fusivities  for  x,  y,  z  (irT2s~-l)  7. 
7,  H:  height  above  the  ground  where  release  occurs  (m)  7. 
7.  U:  wind  velocity  (m/s)  7. 
7.  7. 


7o7.7o7o7o7o7.7o7o7.7o7.7o7o7.7o7o7o7o7.7o7o7.7o7o7o7o7.7o7o7o7o7.7o7o7.7o7o7o7o7.7o7o7o7o7.7o7o7.7o7o7o7o7o7.7o7o7.7o7o7o7o7.7o7o7.7o7o7o7o 


clc ; 
clear; 

format  long; 


7o7.7o7o7o7o7.7o7o7.7o7.7o7o7o7o7o7o7o7o7o7o7.7o7o7o7o7.7o7o7o7o7.7o7o7o7o7o7o7o7o7o7o7o7o7.7o7o7o7o7o7o7o7o7.7o7o7. 


7.  Time  matrix 


7. 


A-l 
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inf 


inf 


inf 

6.29 

inf 

6.52 
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6.29 

inf 
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inf 
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inf 

6.63 

inf 
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6.88 
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inf 

8.4 

inf 

inf 
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6.71 

6.88 

inf 

inf 

inf 

inf 
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inf 

inf 

inf 

inf 

inf 

7 

6.76 

inf 

inf 

inf 

inf 

6.16 

inf 

inf 

inf 

inf 

inf 

inf 

inf 

inf 

9.13 

inf 

inf 

inf 

inf 

inf 

inf 

7.2 

inf 

inf 

inf 

inf 

inf 

7 

inf 

inf 

9.13 

inf 

6.76 

6.16 

inf 

inf 

inf 

inf 

7.18 

10.05 

inf 

inf 

inf 

inf 

inf 

7.43 

7.18 

inf 

inf 

7.94 

5.57 

10.05 

inf 

7.94 

inf 

inf 

inf 

7.43 

5.57 

inf 

inf]  ; 

mmmranfflnmmnnamjimmnnmffl 


’/,  x,  y  coordinates 


y:/.7.y:/:/.y.y:/:/.7.y:/:/.y.7.y:/:/.y.y:/:/.7.y:/:/.y.y:/:/:/.y.y:/;/.7.y:/.7.y.y:/:/.7.y.y:/:/.7.,/:/:/.7.y.y:/.y. 


X  =  [6.275673696 


14.46607868 


33.29676199 


42.99188208 


60.2395703 

71.90038759 

86.36478774 

94.14929655 

98.82869961 

99.741966] ; 


Y  =  [23.01492355 
56.52684103 
96.01718192 
5 . 234534745 
53 . 50444044 
140.5256813 
36.6666158 
86.36143071 
7.663960692 
100.4998932] ; 


y:/.y.y:/:/.y.y:/:/.y.y:/:/.y.y.y:/:/.y.y:/:/.y.y:/:/.y.y:/:/:/.y.y:/:/.y.y.y.y.y.y:/.y.y.y.y:/.y.y.y.y:/:/.y.y:/.y. 


I  Parameters  for  sample  advection-dif fusion  equation  used  '/, 


y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y//.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.y. 


nel  -  10; 


Q  =  1000; 
K_x  =  12; 
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K_y  =  12; 

K_z  =  .2113; 
Z=0; 


H=50 ; 


U= . 25 ; 

V  =  0.0006; 


mmrararara  %mrararammmrarararmmn 


t  =  7*60;  7,  time  since  release 
for  i  =  l:nel 

7,  Computes  the  vapor  concentration  (kg/nT3)  for  location  (x_i,y_i) 
v(i)  =  Q/ (8*pi~  (3/2)  *  (K_x*K_y*K_z)  ~  (1/2)  *t'"  (3/2)  ) 

*exp(- (X(i) -U*t) "2/ (4*K_x*t) -Y (i) ~2/ (4*K_y*t) ) * (2*exp(- (H) "2/ (4*K_z*t) ) ) 
7,  Converts  units  to  mg/nT3 
v(i)  =  v(i)*1000; 

end 

for  i  =  l:nel 

7,  If  the  vapor  concentration  at  location  i  is  below  the  threshold,  a 
7,  reward  value  of  1  is  assigned,  else  a  value  of  zero  is  assigned 
if  v(i)  <  V 
r  (i)  =  1; 

else 

r (i)  =  0; 

end 

end 

v 

r 
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Appendix  B.  Dynamic/Deterministic  Code 


7.  Computing  Dynamic  and  Deterministic  7, 
7,  vapor  concentrations  for  (x,y)  coordinates  7« 
7.  7. 
7,  The  purpose  of  this  MATLAB  code  is  to  compute  the  dynamic  and  l 
7.  deterministic  vapor  concentrations  for  the  network  used  in  7« 
7,  this  thesis.  7« 


7. 


7. 


7,  Author:  Jennifer  R.  Plourde,  AFIT/ENS  G0R-05M  7« 
7,  Date:  1  Feb  05  l 
7.  Last  revised:  1  Feb  05  7« 
7,  References:  Kathirgamanathan,  P.,  McKibbin,  R.  and  % 
7,  R.I.  McLachlan  (2003).  Source  release-rate  estimation  7« 
7,  of  atmospheric  pollution  from  a  non-steady  point  source  % 
7,  -  Part  1:  Source  ate  a  known  location.  Res.  Lett.  Inf.  % 
7.  Math.  Sci.  vol  5,  71-84.  7« 
7.  7. 


7o7.7o7o7o7o7.7o7o7.7o7.7o7o7.7o7o7o7o7.7o7o7.7o7o7o7o7.7o7o7o7o7.7o7o7.7o7o707o7.7o7o7o7o7.7o7o7.7o7o7o7o7o7.707o7.7o7o7o7o7.7o7o7.7o7o7o7o 


nininninnininn  VARIABLE  definitions  mmrammmmra% 


7.  Inputs  required:  7« 
7,  nel:  number  of  nodes  7» 
7,  tm:  (nel  x  nel)  matrix  of  feasible  arc  times  between  nodes  (min)  7# 
7,  X:  (nel  x  1)  vector  of  x-coordinates  °i 
7,  Y:  (nel  x  1)  vector  of  y-coordinates  7# 
7,  Z:  (nel  x  1)  vector  of  z-coordinates  (z  =  0,  for  this  thesis)  % 
7,  vt:  vapor  threshold  value  7« 
7,  Q :  total  mass  release  (kg)  7« 
7,  K_x,  K_y,  K_z:  eddy  dif fusivities  for  x,  y,  z  (m's2s'‘-l)  7« 
7.  H:  height  above  the  ground  where  release  occurs  (m)  7« 
7,  U:  wind  velocity  (m/s)  l 
7,  T:  time  allotted  for  search  7« 
7,  current:  starting  location  7o 
7,  R:  set  of  nodes  in  the  path  l 
7.  reduced:  set  of  nodes  with  reward  of  1  */ 
7,  N:  set  of  nodes  l:nel  7« 
7,  Ni:  set  of  node  adjacent  to  node  i  7« 
7,  tf :  future  time  (t  +  t_i,j)  l 


7o7.7o7o7o7o7.7o7o7.7o7.7o7o7.7o7o7o7o7.7o7o7.7o7o7o7o7o7o7o7o7o7.7o7o7.7o7o7o7o7.7o7o7o7o7.7o7o7.7o7o7o7o7o7.707o7.7o7o7o7o7.7o7o7.7o7o7o7o 


clc ; 


B-l 


clear; 

format  long; 


7:/.y.y////.7.y//:/.y.%y:/:/.y.7:/;/.7.y:/:/.y.y:///:/.y.y;/;/.7.y//:/.y.y:///.%y.y:/;/.7.y:/:/.y.y:/.,/.y.7.,/.y.7.y.7:/. 

I  User  Inputs  ’/, 


y:/.y.y:/:/.y.y:/y/.y.y.y:/:/.y.y:/:/.y.y:/:/.y.y:/:/:/.y.y:/:/.y.y:/:/.y.y:/;/.y.y.y.y.y.y.y.y.y.y.y.y:/.y.y.y:/:/.y.y:/. 


nel  =  10; 

T  =  90; 
vt  =  0.0006; 
current  =  5 ; 
Q  =  1000; 

K_x  =  12; 

K_y  =  12; 

K_z  =  .2113; 
Z=0 ; 

H=50; 

U= . 25 ; 


inf 

6.29 

inf 

6.52 

7.31 

inf 

inf 

inf 

inf 

inf 

6.29 

inf 

6.63 

inf 

6.71 

inf 

inf 

inf 

inf 

inf 

inf 

6.63 

inf 

8.4 

6.88 

7.2 

inf 

inf 

inf 

inf 

6.52 

inf 

8.4 

inf 

inf 

inf 

7 

inf 

inf 

9.13 

7.31 

6.71 

6.88 

inf 

inf 

inf 

6.76 

6.16 

inf 

inf 

inf 

inf 

7.2 

inf 

inf 

inf 

inf 

7.18 

10.05 

inf 

inf 

inf 

inf 

7 

6.76 

inf 

inf 

inf 

inf 

7.43 

inf 

inf 

inf 

inf 

6.16 

7.18 

inf 

inf 

7.94 

5.57 

inf 

inf 

inf 

inf 

inf 

10.05 

inf 

7.94 

inf 

inf 

inf 

inf 

inf 

9.13 

inf 

inf 

7.43 

5.57 

inf 

inf] 

X  =  [6.275673696 
14.46607868 
33.29676199 
42.99188208 
60.2395703 
71.90038759 
86.36478774 
94.14929655 
98.82869961 
99.741966] ; 

Y  =  [23.01492355 
56.52684103 
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96.01718192 


5 . 234534745 

53 . 50444044 

140.5256813 

36.6666158 

86.36143071 

7.663960692 

100.4998932] ; 

reduced  =  []  ; 

R  =  [current] ; 


for  i  =  l:nel 
N(i)  =  1; 

end 

7  initial  starting  time 
t  =  7*60;  '/  in  seconds 

"/Initialization  step 
for  i  =  current 

v(i)  =  Q/ (8*pi~ (3/2) * (K_x*K_y*K_z) ~ (1/2) *t~ (3/2) ) 

*exp(-  (X ( i )  -U*t)  ~2/  (4*K_x*t)  -Y (i)  ~2/  (4*K_y*t)  )  *  (2*exp(-  (H)  '‘2/  (4*K_z*t)  )  ) 
v(i)  =  v(i)*1000;  "/converts  to  mg 
if  v(i)  <  vt  &  v(i)  >  0 
r(i)  =  1; 

else 

r (i)  =  0; 

end 

end 


*/  If  initial  v  <  V,  current  is  subtracted  from  set  of  nodes  to  visit 
7,  Otherwise,  it  is  still  eligible  to  be  searched 

N(i)  =  N(i)  -  1; 


7o7omr/.7oEnd  initialization mmrammramm 


mmmrararammmrarammmmnm 


7,  Step  1 


7. 


"/"/"/“/7."/7o7.“/7o7o“/*/7»7o7.7o7.7o7«7o7.7o7o7.7."/"/*/“/7o7«7«7.7o7o7«7.*/"/"/“/7o7o7o7o7.7o7o7o 


b  =  0; 

while  t  <  90*60 
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N; 

for  f  =  l:nel 

if  tm( current,  f)  ==  inf 
Ni (f )  =  0; 

else 

Ni (f )  =  1; 

end 

end 

Ni; 

t  =  t/60;  '/return  time  units  to  minutes  because  of  tm 

'/Calculate  vapor  concentration  for  each  node  j 
for  j  =  l:nel 

tf  =  t  +  tm (current,  j); 

tf  =  tf*60;  '/converts  to  seconds 

v(j)  =  Q/ (8*pi~ (3/2) * (K_x*K_y*K_z) ~ (1/2) *tf ~ (3/2) ) 

*exp(- (X( j ) -U*tf ) ~2/ (4*K_x*tf ) -Y ( j ) ~2/ (4*K_y*tf ) ) * (2*exp(- ( (H) ~2) / (4*K_z*tf ) ) ) ; 
v(j)  =  v(j)*1000;  7,  converts  to  mg 

end 

for  j  =  l:nel 
if  v(j)  <  0 

v(j)  =  -v(j); 

end 

end 

false  =  0; 
used  =  false; 

7#  Allow  return  to  previously  visited  node  if  there  is  no  other  option 
if  max(N  +  Ni)  <=  1 
ntmp  =  N; 

N  =  N  +  Ni; 
used  =  1; 

end 

7,  If  a  node  has  been  visited,  set  v  <  0  so  it  won’t  be  considered 
for  j  =  l:nel 

if  Ni ( j )  ==  1 
v(j)  =  v(j); 
elseif  Ni(j) 

v(j)  =  -v(j); 

end 
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end 


if  b  >  2 

temp  =  R(b  -  1) ; 

if  v(temp)  >  0 

v(temp)  =  -v(temp) ; 

end 

end 


if  used 

N  =  ntmp; 

end 

v 

°/0  If  v  <  V,  set  r  =  1  (MOPP  reduced)  o.w.  r  =  0 
for  j  =  l:nel 

if  v(j)  <  vt  &  v(j)  >=  0 
r(j)  =  1; 

else 

r(j)  =  0; 

end 


end 


7o7oy//.7.7o7.%7.7.7o0/.%%7.yoyo7.%%yo7o%%7//.7o%%%7.7o7.%%7.yoy.7.%7.yoyo%%%0/.7o%% 


"/,  Step  2 


7. 


y.y.y.y.y.y.y.y.y.y.%y.y.y.y.y.y.y.%y.y.y.y.y.y.y.%y.y.y.y.y.y.y.y.y.y.y.y.%y.y.y.y.y.y.y.y.y.y. 


if  r  ==  0 

for  j  =  l:nel 
if  v(j)  <  0 

ty(j)  =  inf; 

else 

ty(j)  =  tm ( current ,j ) ; 


end 

M  =  min(ty)  ; 

end 

end 

if  sum(r)  >  0 

for  j  =  l:nel 

if  r(j)  ==  1 

M(j)  =  v(j ) ; 
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else 

M(j)  =  0; 

end 

end 

end 

P  =  max(M) ; 
k  =  0; 

*/,  Record  j  as  the  next  node 
for  j  =  l:nel 

if  v ( j )  ==  P 
next  =  j  ; 

elseif  tm(current,  j)  ==  P 
next  =  j ; 

end 

end 

for  j  =  l:nel 

if  next  ==  j 

if  r(j)  ==  1 
k  =  j; 

end 

end 

end 

if  k  “=  0 

reduced  =  [reduced,  k] ; 

end 

°/0  Remove  j  from  further  consideration 
for  j  =  l:nel 

if  j  ==  next 

7.  r  ( j  )  =  1; 

N( j )  =  N( j )  -  1; 

end 

end 

"/update  time  step 

t  =  t  +  tm(current,  next)  ;  */  minutes 
t  =  t*60;  ‘/converts  to  seconds 

if  r(next)  ==  1 
N(next)  =  -1; 

end 
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°/0Set  next  node  as  current  and  return  to  step  1 


current  =  next ; 
if  N  ==  0 

for  j  =l:nel 

if  r  ( j  )  ==  0 

N(j)  =  1; 

else 

N ( j )  =  0; 

end 

end 

end 

R  =  [R,  next] ; 
b  =  length (R) ; 
if  N  <  0 
break 

end 

ed  =  length (R) ; 

edr  =  length (reduced) ; 


°/0  Step  3 


•/. 


my.y.mmmmmmmmy.mm,mmy.m%m 


if  t  >  90*60 


over  =  length (R)  -  1; 

R(over) ; 

time  =  t  -  tm(R(over) ,R(length(R)))*60; 

R(ed)  =  R(ed)  -  R(ed) ; 

reduced(edr)  =  reduced(edr)  -reduced(edr) ; 
for  j  =  l:nel 

v(j)  =  Q/ (8*pi~ (3/2) * (K_x*K_y*K_z) ~ (1/2) *t~ (3/2) ) 

*exp(- (X( j ) -U*t) ~2/ (4*K_x*t) -Y ( j ) ~2/ (4*K_y*t) ) * (2*exp(- ( (H) ~2) / (4*K_z*t) ) ) ; 
v(j)  =  v(j)*1000;  °/0mg 
for  j  =  l:nel 

if  v(j)  <  vt 
r(j)  =  1; 

else 

r  (j  )  =  0; 

end 

end 
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end 


end 

end 

path  =  R; 
path 

edr  =  edr  -1; 

disp([’The  number  of  areas  in  reduced  MOPP  is  ’  int2str (edr)] ) 
disp(  [’These  areas  are  ’  int2str (reduced)] ) 

disp([’The  total  time  of  the  search  is  ’  num2str (time/60)  ’  minutes.’]) 
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Appendix  C.  Dynamic/Stochastic  Code 


7:/.y.y//:/.7.y:/:/.y.7.y:/:/.y.7:/:/.7.y:/:/.y.7:/:/:/.y.7:/:/.7.y:/:/.y.y:/;/.7.y.y:/;/.7.y:/.7.y.y:/:/:/.7.,/:/.7.y.7:/:/.7.,/:/.7.y.y:/:/. 


7.  Computing  Dynamic  and  Stochastic  7# 
7,  vapor  concentrations  for  (x,y)  coordinates  7« 
7.  7. 
7,  The  purpose  of  this  MATLAB  code  is  to  compute  the  dynamic  and  7» 
7.  stochastic  vapor  concentrations  for  the  network  used  in  7« 
7,  this  thesis.  7# 


7. 


7. 


7,  Author:  Jennifer  R.  Plourde,  AFIT/ENS  G0R-05M  7« 
7,  Date:  1  Feb  05  l 
7.  Last  revised:  1  Feb  05  7« 
7,  References:  Kathirgamanathan,  P.,  McKibbin,  R.  and  % 
7,  R.I.  McLachlan  (2003).  Source  release-rate  estimation  7« 
7,  of  atmospheric  pollution  from  a  non-steady  point  source  % 
7,  -  Part  1:  Source  ate  a  known  location.  Res.  Lett.  Inf.  % 
7.  Math.  Sci.  vol  5,  71-84.  7« 


7o7.7o7o707o7.7o7o7.7o7o7o7o7.7o7o7o7o7o707o7.7o7o7o7o7o7o7o707o7.7o7o7o7o7o7o7o7«7o7o7o7o7.7o7o7o7o7o7o7o7o7.707o7o7o7o7o7o7.707o7o7o7o707o 


7.7.7o7o707o7.7o7o7o7o7.7o7o7.7o7o7o7o7.7o7o7.  variable  definitions  llliniUiniUinnillll 


7,  Inputs  required:  7« 
7,  nel:  number  of  nodes  7« 
7.  tm:  (nel  x  nel)  matrix  of  feasible  arc  times  between  nodes  (min)  % 
7,  X:  (nel  x  1)  vector  of  x-coordinates  7« 
7,  Y:  (nel  x  1)  vector  of  y-coordinates  7# 
7,  Z:  (nel  x  1)  vector  of  z-coordinates  (z  =  0,  for  this  thesis)  % 
7,  vt:  vapor  threshold  value  7« 
7,  Q :  total  mass  release  (kg)  l 
7,  K_x,  K_y,  K_z:  eddy  dif fusivities  for  x,  y,  z  (nT2s~-l)  7« 
7,  H:  height  above  the  ground  where  release  occurs  (m)  7« 
7,  U:  wind  velocity  (m/s)  7o 
7.  T:  time  allotted  for  search  / 
7,  current:  starting  location  7« 
7,  R:  set  of  nodes  in  the  path  7« 
7,  reduced:  set  of  nodes  with  reward  of  1  °( 
7,  N:  set  of  nodes  l:nel  7» 
7.  Ni:  set  of  node  adjacent  to  node  i  7« 
7,  tf :  future  time  (t  +  t_i,j)  7« 
7,  mu:  (lx  nel)  vector  of  expected  value  of  7« 
7,  vapor  concentration  for  each  node  7« 


7o7.7o7o7o7o7.7o7o7.7o7.7o7o7.7o7o7o7o7.7o7o7.7o7o7o7o7o7o7o7o7o7.7o7o7.7o7o7o7o7.7o7o7o7o7.7o7o7.7o7o7o707o7.707o7.7o7o7o7o7.7o7o7.7o7o7o7o 


clc ; 
clear; 
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format  long; 


nel  =  10; 

T  =  90; 
vt  =  0.0006; 
current  =  5 ; 

Q  =  1000; 

K_x  =  12; 

K_y  =  12; 

K_z  =  .2113; 

Z=0 ; 

H=50; 

U= . 25 ; 

mu  =  [  .00015  .00065  .000756  .000635  .00057  .0001  .00023  .00045  .0011  .00013]  j’/.EV’s 


tm  =[inf 

6.29 

inf 

6.52 

7.31 

inf 

inf 

inf 

inf 

inf 

6.29 

inf 

6.63 

inf 

6.71 

inf 

inf 

inf 

inf 

inf 

inf 

6.63 

inf 

8.4 

6.88 

7.2 

inf 

inf 

inf 

inf 

6.52 

inf 

8.4 

inf 

inf 

inf 

7 

inf 

inf 

9.13 

7.31 

6.71 

6.88 

inf 

inf 

inf 

6.76 

6.16 

inf 

inf 

inf 

inf 

7.2 

inf 

inf 

inf 

inf 

7.18 

10.05 

inf 

inf 

inf 

inf 

7 

6.76 

inf 

inf 

inf 

inf 

7.43 

inf 

inf 

inf 

inf 

6.16 

7.18 

inf 

inf 

7.94 

5.57 

inf 

inf 

inf 

inf 

inf 

10.05 

inf 

7.94 

inf 

inf 

inf 

inf 

inf 

9.13 

inf 

inf 

7.43 

5.57 

inf 

inf]  ; 

X  =  [6.275673696 
14.46607868 
33.29676199 
42.99188208 
60.2395703 
71.90038759 
86.36478774 
94.14929655 
98.82869961 
99.741966] ; 

Y  =  [23.01492355 
56.52684103 
96.01718192 
5 . 234534745 
53 . 50444044 
140.5256813 


C-2 


36.6666158 


86.36143071 
7.663960692 
100.4998932] ; 

reduced  =  []  ; 

R  =  [current] ; 

for  i  =  l:nel 
N(i)  =  1; 

end 

°/0Time  since  release 
t  =  7*60  ;°/0seconds 

rand( ’ state ’ ,  135); 


mmmmrammmmmmmmraram 


°/0  Initialization 


7. 


mmraram%mrarararammmraray.ra 


for  i  =  current 

V(i)  =  expcdf(vt,  mu(i)); 
v(i)  =  exprnd(mu(i) ) ; 
if  v(i)  <  vt  &  v(i)  >  0 
r  (i)  =  1; 

reduced  =  [reduced,  i] ; 


else 


r (i)  =  0; 

end 


end 


N(i)  =  N(i)  -  1; 

7o7.7.y.7.7.7.7oEnd  Initialization%%%%mm%%%%%%%mm% 

mmmrararammmmmrammranm 


7,  Step  1 


7. 


mmrarararammmrararammmnm 


b  =  0; 

while  t  <  90*60 
for  f  =  l:nel 

if  tm ( current ,  f)  ==  inf 
Ni (f )  =  0; 

else 

Ni (f )  =  1; 
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end 


end 

t  =  t/60;  yoreturn  time  units  to  minutes  because  of  tm 

for  j  =  l:nel 

tf  =  t  +  tm (current,  j); 
tf  =  tf*60; 
if  Ni ( j )  ==  1; 

mum(j)  =  mu( j ) * . 01*tf ; 

V ( j )  =  expcdf(vt,  mum(j)); 

else 

V(j)  =  0; 

end 

end 

V 

for  j  =  l:nel 
if  V(j)  <  0 

V(j)  =  -V(j); 

end 

end 


false  =  0; 
used  =  false; 

°/0  Allow  return  to  previously  visited  node  if  there  is  no  other  option 
if  max(N  +  Ni)  <=  1 
ntmp  =  N; 

N  =  N  +  Ni; 
used  =  1; 

end 


rararammmrararammmmrararan 


°/0  Step  2 


°/o 


rara%mrarararammmrararayx/.mm% 


l  If  a  node  has  been  visited,  set  v  <  0  so  it  won’t  be  considered 
for  j  =  l:nel 

if  Ni ( j )  ==  1 
V(j)  =  V( j )  ; 
elseif  Ni ( j ) 

V(j)  =  — V ( j ) ; 

end 


end 


if  b  >  2 

temp  =  R(b  -  1) ; 
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if  V(temp)  >  0 

V(terap)  =  -V(temp); 

end 

end 

if  used 

N  =  ntmp; 

end 

P  =  max(V) ; 
mts  =  []  ; 

for  j  =  l:nel 

if  Ni ( j )  ==  1; 
if  V(j)  ==  P 

mt  =  tm( current , j) ; 

if  mt  >0 

mts  =  [mts,  mt] ; 

if  tm (current,  j) 
next  =  j  ; 

end 

end 

end 

end 

end 

for  i  =  l:nel 

v(next)  =  exprnd(mu(next) ) ; 

end 

for  j  =  next 

if  v(j)  <  vt  &  v(j)  >=  0 
r(j)  =  1; 

else 

r(j)  =  0; 

end 

end 

if  r(next)  ==  1 
k  =  next ; 

reduced  =  [reduced,  k] ; 


==  min (mts) 
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end 


t  =  t  +  tm( current,  next)  "/minutes 

t  =  t*60; 

next 

if  r(next)  ==  1 
N(next)  =  -1; 

end 

current  =  next ; 
if  N  ==  0 

for  j  =l:nel 

if  r(j)  ==  0 
N(j)  =  1; 

else 

N(j)  =  0; 

end 

end 

end 

R  =  [R,  next] ; 
b  =  length (R)  ; 
v 

if  N  <  0 
break 

end 

ed  =  length (R) ; 

edr  =  length (reduced) ; 

if  t  >  90*60 

R(length(R) ) ; 

over  =  length (R)  -  1; 

time  =  t  -  tm(R(over) ,R(length(R) ) ) *60 ; 

R(ed)  =  R(ed)  -  R(ed) ; 

reduced(edr)  =  reduced(edr)  -  reduced(edr) ; 

end 

end 

path  =  R; 
path 

S  =  sum(r) ; 

disp([’The  number  of  areas  in  reduced  MOPP  is  ’  int2str(S)]) 
disp(  [’These  areas  are  ’  int2str (reduced)] ) 

disp([’The  total  time  of  the  search  is  ’  num2str (time/60)  ’  minutes.’]) 
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Appendix  D.  Static  Cases  Enumeration  Code 


7:/.y.y//:/.7.y:/:/.y.7.y:/:/.y.7:/;/.7.y:/:/.y.7:/:/:/.y.y:/;/.7.y:/:/.y.y:/;/.7.y.y:/;/.7.y:/.7."/.y:/:/:/.7.,/:/.7.y.7:/:/.7.,/.7.7.y.y:/:/. 


7.  Enumeration  of  all  feasible  paths  in  static/deterministic  and  / 


%  static/stochastic  case  l 

7.  7. 
7,  The  purpose  of  this  code  is  to  find  the  optimal  shortest  % 
7.  path  of  maximum  reward  for  the  static/deterministic  and  7« 
7,  static/stochastic  case  of  this  thesis.  7« 
7.  7. 
7.  Author:  Jennifer  R.  Plourde,  AFIT/ENS  G0R-05M  7« 
7.  Date:  10  Feb  05  7. 
7.  Last  revised:  11  Feb  05  7« 
7,  References:  None  7» 


7o7.7o7o7o7o7.7o7o7.7o7.7o7o7.7o7o7o7o7.7o7o7.7o7o7o7o7.7o7o7o7o7.7o7o7.7o7o707o7.7o7o7o7o7.7o7o7.7o7o7o7o7o7.707o7.7o7o707o7.7o7o7.7o7o7o7o 


7o7.7o7o7o7o7o7o7o7o7o7o7o7o7o7o7«7o7o7.7o7«7.  variable  definitions  llinillllllUllUllllUl 


7.  Inputs  required:  7« 

7,  tm:  (nel  x  nel)  matrix  of  feasible  arc  times  between  nodes  (min)  7« 


7o7.7o7o7o7o7.7o7o7.7o7.7o7o7.7o7o7o7o7.7o7o7.7o7o707o7.7o7o7o7o7.7o7o7.7o7o707o7.7o7o7o7o7.7o7o7.7o7o7o7o7o7.707o7.7o7o7o7o7.7o7o7.7o7o7o7o 


clc ; 
clear; 
tO  =  clock; 


[inf 

6.29 

inf 

6.52 

7.31 

inf 

inf 

inf 

inf 

inf 

6.29 

inf 

6.63 

inf 

6.71 

inf 

inf 

inf 

inf 

inf 

inf 

6.63 

inf 

8.4 

6.88 

7.2 

inf 

inf 

inf 

inf 

6.52 

inf 

8.4 

inf 

inf 

inf 

7 

inf 

inf 

9.13 

7.31 

6.71 

6.88 

inf 

inf 

inf 

6.76 

6.16 

inf 

inf 

inf 

inf 

7.2 

inf 

inf 

inf 

inf 

7.18 

10.05 

inf 

inf 

inf 

inf 

7 

6.76 

inf 

inf 

inf 

inf 

7.43 

inf 

inf 

inf 

inf 

6.16 

7.18 

inf 

inf 

7.94 

5.57 

inf 

inf 

inf 

inf 

inf 

10.05 

inf 

7.94 

inf 

inf 

inf 

inf 

inf 

9.13 

inf 

inf 

7.43 

5.57 

inf 

inf] 

[jpath, jtime]  =  jenndykstra(tm) 
tl  =  clock; 
etime (tl , tO) 
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mmmxmmxmm  Main  code  mmrara%%mmmram%ram 

7,  Inputs  required:  7« 

7,  tm:  (nel  x  nel)  matrix  of  feasible  arc  times  between  nodes  (min)  *( 

7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7o7.7o7o7o7o7o7o7o7o7o7.7o7o7o7o7o7o7o7p7o7o7o7o7p7o7o7o7o7o7o7»7o7.7o7o7o7o7o7o7o7p7o7o7o7o7o7o7o 


function  [ j path, j time]  =  jenndykstra(tm) 

^building  vectors 

leftv  =  ones(l , size(tm, 1) ) ; 

7pStart  at  node  5 

currv  =  [zeros (1,  4),  1,  zeros (1,5)]; 

7pcurrv  =  [l,zeros(l,size(tm,  1)-1)]  ; 
prevv  =  zeros(l,size(tm,l)) ; 
maxtime  =  90; 
numbervector  =  []  ; 
for  i  =  l:size(tm,l) 

numbervector  =  [numbervector , i] ; 

end; 

streamv  =  [] ;  %  current  path 

tims  =0;  l  accumulated  time 

7odet ermine  maximum  size  of  path 

maxsize  =  ceil(maxtime/min(min(tm)))+20; 

^create  tree  and  associated  times 

[masters ,mastert]  =  steptree3(tm,  leftv,  currv,  prevv,  streamv,  tims, 
numbervector,  maxsize); 

7#choose  shortest  time  from  the  list 
shortest  =  min(mastert) ; 
found  =  0; 
i  =  1; 

7.  put  rewards  here 

^identify  the  path  associated  with  the  shortest  time 
while  found  <  1 

if  mastert(i,l)  ==  shortest 
jpath  =  masters (i, :) ; 
jtime  =  mastert (i , : ) ; 
found  =  1; 

end; 

i  =  i  +  1; 


maxtime , 
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end; 


i  =  1; 

while  ((jpath(l,i)  ~=  0)  &  (i  <  maxsize)) 
i  =  i  +  1; 

end; 

jpath  =  jpath(l , 1 : i-1) ; 


7«°/0  this  is  the  subroutine  that  is  called  recursively 


function  [masterss,  mastertt]  =  steptree3(timem,  leftv,  currv,  prevv,  streamy,  tims, 

numbervector ,  maxsize) 


masterss  =  []  ; 
mastertt  =  []  ; 

'/update  step 
if  leftv*currv’  ==  1 

leftv  =  leftv-currv; 

end; 

tempstreamv  =  streamv; 

streamv  =  [streamv, currv*numbervector ’] ;  7  where  we’re  at 

leftv  =  logical (leftv) ; 
currv  =  logical (currv) ; 
prevv  =  logical (prevv) ; 

if  sum(prevv)  ~=  0 
temptims  =  tims; 

tims  =  tims  +  t imem (prevv, currv) ; 

end; 

if  ( (sum(leftv)  ==  0)  &  (  tims  <=  maxtime)) 

7  if  ( (sum (leftv)  ==  0)  I  (  tims  >=  maxtime)) 
masterss  =  [streamv, zeros(l,maxsize-size(streamv, 2))] ; 
mastertt  =  tims; 
elseif  (  tims  >  maxtime) 

masterss  =  [tempstreamv,zeros (1 ,maxsize-size (tempstreamv, 2))] ; 
mastertt  =  temptims; 


else 

tempstring  =  timem( :, currv) ; 
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maxt ime , 


run  =  0; 


for  i  =  1 : size (timem, 1) 

if  ( (tempstring (i , 1)  <=  maxtime)  &  (leftv(l,i)  ==  1)) 
run  =  1; 

nextstep  =  zeros (1 , size (timem, 1) ) ; 
nextstep(l , i)  =  1; 

[temps,  tempt]  =  steptree3 (timem,  leftv,  nextstep,  currv,  streamv,  tims,  maxtime, 
numbervector ,  maxsize) ; 
if  size (tempt , 1)  >=1 

masterss  =  [masterss ; temps] ; 
mastertt  =  [mast ertt ; tempt] ; 

end; 

end; 

end; 

if  run  ==  0  “/*  called  when  stuck 
for  i  =  1 : size(timem, 1) 

if  tempstring (i , 1)  <=  maxtime 

nextstep  =  zeros (1, size (timem, 1)) ; 
nextstep (1 , i)  =  1; 

[temps,  tempt]  =  steptree3 (timem,  leftv,  nextstep,  currv,  streamv,  tims,  maxtime, 
numbervector,  maxsize); 
if  size (tempt , 1)  >=1 

masterss  =  [masterss ; temps] ; 
mastertt  =  [mastertt ; tempt] ; 

end; 

end; 

end; 

end; 

end; 
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